Analyses for paper

# colors for gender
sex_colors =  c(men = "#30D5C8",women = "purple",
                `women-men` = "orange", avg = "blue4")
scale_col_sex = list(
  scale_fill_manual(values = sex_colors), 
  scale_color_manual(values = sex_colors)  
)

sex_colors.b =  c(men = "#30D5C8",women = "purple",
                `Sex difference` = "orange", total = "blue4")
scale_col_sex.b = list(
  scale_fill_manual(values = sex_colors.b), 
  scale_color_manual(values = sex_colors.b)  
)

sex_colors2 =  c("#30D5C8","purple")
scale_col_sex2 = list(
  scale_fill_manual(values = sex_colors2), 
  scale_color_manual(values = sex_colors2)  
)
# colors for drug
drug_colors =  c("blue","grey25","orange")
scale_col_drug = list(
  scale_fill_manual(values = drug_colors), 
  scale_color_manual(values = drug_colors)  
)
# colors for stress condition
stress_colors =  c("#1b9e77","#d95f02")
scale_col_stress = list(
  scale_fill_manual(values = stress_colors), 
  scale_color_manual(values = stress_colors)  
)
gg_no_y_axis = 
  theme(axis.text.y = element_blank(),
        axis.ticks.length.y = unit(0,"mm"),
        panel.grid.major.y = element_blank()) 
gg_no_x_axis = 
  theme(axis.text.x = element_blank(),
        axis.ticks.length.x = unit(0,"mm"),
        panel.grid.major.x = element_blank()) 

negaffect_items = c("distressed","anxious","vulnerable","irritable") 
posaffect_items = c("good", "happy", "confident", "safe", "relaxed")
embselfcn_items = c("embarassed","selfconscious")
sideffect_items = c("dizzy","dry_mouth", "blunted", "nauseous","warm_face","not_myself") 

1 Prepare data

We load the data. The data file was created on NA.

load("data/my_data_Questionnaire.Rdata")
my_data.Q = selector(my_data)
my_data.Q[item_name == "dull", item_name := "blunted"]

load("data/my_data_cortisol.Rdata")
my_data.cortisol = 
  my_data %>% 
  .[, result := result*100*0.0001] %>% 
  .[, log.result := log(result)] %>% 
  .[participant %in% unique(my_data$participant)]

load("data/my_data_SelfAdmin.Rdata")
my_data = selector(my_data)

stage.names = c("BL 1", "BL 2", "BL 3", "state induction", "dose 1" , 
           "reminder 1", "self- admin", "reminder 2", "DM task",
           "cognitive task", "dose 2") %>% gsub(" ","\n",.)
my_data %>% 
  .[, state.condition := ifelse(Stress == "stress", "stress", "control")] %>% 
  .[, Sex := ifelse(gender == "m","men","women")] 
my_data.Q %>% 
  .[, state.condition := ifelse(Stress == "stress", "stress", "control")] %>% 
  .[, Sex := ifelse(gender == "m","men","women")] %>% 
  .[, induced.state := ifelse(Stress == "stress" & stage > 3, "stress", "control")] %>% 
  .[, Stage := factor(stage.names[stage], levels = stage.names)]

Calculate first time at target for self admin:

my_data %>% 
  .[, first.time.target := 1e5 + 0.1] %>% 
  .[, `:=`(first.time.target = 
          get_ftt(effect_obtained,target, time, bias),
          mean.delta_effect_obtained = mean(abs(effect_obtained-target))),
        by = .(participant,session)]

Currently 63 participants (32 women) with data from 4 opioid sessions are included in the analysis.

1.1 Exclusions based on ratings

We exclude people who reported neither an increases in negative emotions nor a decrease in positive emotions after the stress condition.

negaffect_items = c("distressed","anxious","vulnerable","embarassed","irritable") ## add stressed

posaffect_items = c("good", "happy", "confident", "safe", "relaxed")

# maximum change in any of the negative emotion items
my_data.negaffect = 
  my_data.Q %>% 
  .[stage %in% c(3,4) & item_name %in% negaffect_items & state.condition == "stress" & Drug == "oxycodone", 
    .(participant, item_name, response, stage, Sex)] %>% 
  dcast( participant + item_name + Sex ~ stage, value.var = "response") %>% 
  .[, delta := `4` - `3`] %>% 
  .[, .(max.neg_incr = max(delta)), by = .(participant, Sex)]

# maximum change in any of the positive emotion items
my_data.posaffect = 
  my_data.Q %>% 
  .[stage %in% c(3,4) & item_name %in% posaffect_items & state.condition == "stress" & Drug == "oxycodone", 
    .(participant, item_name, response, stage, Sex)] %>% 
  dcast( participant + item_name + Sex ~ stage, value.var = "response") %>% 
  .[, delta := `4` - `3`] %>% 
  .[, .(max.pos_decr = min(delta)), by = .(participant, Sex)]


my_data.posng_change = 
  merge(my_data.negaffect,
        my_data.posaffect,
        by = c("participant","Sex")) %>% 
  setkeyv(c("participant")) %>% 
  .[, stress.filter := ifelse(max.neg_incr <= 0 & max.pos_decr >= 0,1,0)]

my_data.posng_change %>% 
  ggplot(aes(x = max.neg_incr, y = max.pos_decr, color = Sex, label = participant)) + 
  geom_text() + 
  theme(legend.position = "right") +
  scale_col_sex +
  ylab("mac increase negative emotion items") +
  xlab("mac decrease positive emotion items")

Figure 1.1: Max changes in self-reported negative and positive items from pre to post stress induction.

0 participants had to be excluded due to insensitivity to the stress induction.

my_data %>% 
  .[my_data.posng_change, stress.filter := stress.filter] 
my_data = my_data[stress.filter == 0]

my_data.Q_keys = key(my_data.Q)
my_data.Q %>% 
  setkeyv(c("participant")) %>% 
  .[my_data.posng_change, stress.filter := stress.filter] %>%  
  setkeyv(my_data.Q_keys)
my_data.Q = my_data.Q[stress.filter == 0]

rm(my_data.posaffect,my_data.posng_change, my_data.negaffect)

Some overview over participants:

• number participants with self admin data: 63
• number of missing sessions from self admin data: 0
• number participants with questionnaire data: 63
• number of missing sessions from questionnaire data: 1

1.2 Following instructions in the self administration task

We check if participants actually work to obtain the stated desired effect. Participants who end at more than 50 points away from the their declared goal are considered to not have followed the instructions. For these participants the questionnaire data are removed from the analysis.

my_data %>% 
  .[, delta_target := abs(tail(effect_obtained,1) - target[1]), by = .(participant, session)] %>% 
  .[, cheater := delta_target > 50]

tmp = 
  my_data[, .(participant,state.condition,session,Drug, cheater)] %>% 
  unique()

cheaters = 
  tmp %>% 
  .[cheater == TRUE] %>% 
  unique() %>% 
  .[, session := factor(session)] 

cheaters  %>% 
  copy() %>% 
  setnames(c("cheater"),c("Ignored instructions")) %>% 
  .[Drug == "oxycodone"] %>% 
  kable() %>% 
  add_footnote(paste("In",tmp[Drug == "placebo" & cheater == TRUE] %>% nrow(),"of the placebo sessions a participant ignored the instructions.")) %>% 
  kable_styling(full_width = FALSE)

participant	state.condition	session	Drug	Ignored instructions
14	control	3	oxycodone	TRUE
14	stress	4	oxycodone	TRUE
15	control	4	oxycodone	TRUE
33	control	3	oxycodone	TRUE
33	stress	4	oxycodone	TRUE
63	control	3	oxycodone	TRUE
63	stress	4	oxycodone	TRUE
77	stress	1	oxycodone	TRUE
77	control	2	oxycodone	TRUE
97	control	1	oxycodone	TRUE
120	stress	3	oxycodone	TRUE
a In 14 of the placebo sessions a participant ignored the instructions.

In 11 of the () oxycodone sessions the obtained effect was more than 50 points from the target away. Altogether 7 participants contributed such sessions. In 14 of the placebo sessions a participant ignored the instructions.”

tmp %>% 
  .[, .(n = sum(!cheater), p = mean(!cheater)), 
    by = .(Drug, state.condition)] %>%
  .[, N := paste0(n, " (", round(p*100),"%)")] %>% 
  dcast(Drug ~ state.condition, value.var = "N") %>% 
  kable(caption = "Sample size for analyses of drug wanting after removal of participants who did not follow instructions") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic() %>%  
  add_footnote("The number in parentheses indicates the percentage of participant who followed instructions.", notation = "none")

Table 1.1: Sample size for analyses of drug wanting after removal of participants who did not follow instructions

Drug	control	stress
oxycodone	57 (90%)	58 (92%)
placebo	55 (87%)	57 (90%)
The number in parentheses indicates the percentage of participant who followed instructions.

2 Analysis plan

2.1 Confirmatory analyses

Due to the within-participant manipulation of the independent variables (drug * state) and the non-normal distributions of several of the outcome variables we will use hierarchical Bayesian regression analyses.

In these analyses, random effects will account for repeated measurement of participants (and time whenever measures are repeated within participant). To account for the non-normal distribution of outcome variables (preliminary analysis of pilot data showed some heteroscedasticity of residual from a gaussian / normal model) we will use an item response theory approach to model the data when appropriate.

See the sections Linear regression and Linear regression without sex effects for an analyses of predictions and residuals of linear models that shows that a linear models are not appropriate for rating data or behavioral data from the self administration task. Note that all outcomes are measured on continuations scales, and we use ordinal logistic regressions merely as a flexible modelling approach that allows us to accurately capture patterns in the data. Because the analysed data is generally continuous, it is possible and permissible to calculate expected responses and derived contrasts from the results of an ordinal logistic regression model.

As a preparation for using and ordinal model, we coarsen the data by binning the 0-100 VAS scale into 20 bins of 5 points and the data of the self-administration task, which has a scales of 0-125, into 25 bins or response levels. The motivation for this bin sizes is that it allows to represent questionnaire responses and behavior relatively fine-grained while being coarse enough to allow for relatively easy estimation of threshold parameters of ordinal logistic models.

We use the more flexible adjacent category ordered logistic model because this allows to relax the proportional odds assumptions if needed by estimating category-specific effects. The priors for the thresholds for this model are set as normal distributions with the expected threshold given no experimental effects as mean and2 as the standard deviation. For examples, given four ordered categories with frequencies 10, 20, 40, 30, the expected threshold between categories 1 and 2 is log((10/70)/(20/70)) = -0.3 and the expected threshold between categories 4 and 4 is log((40/70)/(30/70)) = 0.125. Additionally, we use weakly informative priors of N(0,2) for the standard deviation of random effects and N(0,2). We also use weakly informative priors for fix effects of interest. See the section Prior sensitivity analysis in the Main Results for a prior sensitivity analysis for the main outcome.

All models are estimated with the R package brms, which uses the Hamiltonian Markov Chain Monte Carlo Sampler implemented in the programming language Stan to obtain posterior distribution of model parameters. Sampling is generally done with 4 chains, each of which uses 1000 samples for warm up and 1000 post-warmup samples. Warmup samples are discarded and all reported results rely on post-warmup samples. We verify model convergence by checking that Rhat values are below 1.001 and the the effective sample size is above 500 (it is usually above 2000). If 1000 post-warmup samples were not sufficient to achieve 500 effective samples, we re-ran the model with 2000 post-warmup samples.

We report all results by calculating effects (means and 95% credible intervals) on the original scale. We describe effects whose credible intervals do not include 0 as “present” and report in addition to means and credible intervals the posterior probability that an effect is larger than 0 (\(P_{E>0}\)). For ordered logistic regressions we calculate effects as follows: (1) We calculate the average response \(\mu_l\) for each category \(l\) of the coarsened variable from the raw data. (2) Using the posterior distribution from the fitted model, we calculate for each condition \(k\) the level-probabilities \(P_{l,k}\), i.e. the probability with which a response in each category is expected. ¹ (3) We calculate the expected response \(R_k\) for each condition as the weighted average of \(mu\) with the level-probabilities \(P_{l,k}\) as weights. (4) We calculate effects of interest as contrasts between the expected responses in different conditions. For instance, stress effects are generally calculated as changes in stress response from pre to post state induction.

This regression analysis also allows adjusting for covariates such as session, drug and state induction order, gender, age, BMI, body weight, time of day, experimenter/medical staff identity, drug side effects.

The our analyses adjust for session number, whereby we model the categorical sessions variable hierarchically, in order to regularize session effects. Given that preliminary analyses showed strong sex effects, we also modeled main and interaction effects for sex. For these analysis we generally calculated sex-specific and over-all (averaged over sex) effects as well as sex differences. However, consistent with the pre-registration, we also estimate a model without an interaction with sex for the main outcome from the self administration task. Where relevant, we model repeated measurement of the same item with random effects.

An initial check of the data also revealed substantial effect heterogeneity, i.e. there was great variation in the response to stress and drug manipulations. We took this into account by explicitly modeling it as random treatment effects. Indeed, model comparison of models with and without random treatment effects showed that the latter has a better out of sample predictive value.

Main predictors: Drug condition, stress condition (Note that some of the hypotheses listed below are tested within one of the drug or stress conditions only).

Please see the supplementary materials for complete specifications of all analyses and code that can be used to re-produce the results.

3 Manipulation check

3.1 Stress

3.1.1 Linear regression

my_data.Q.stressed = select_data("stressed")
my_stages = c(1:5)

my_data.Q.stressed %>% plot_stages()

Figure 3.1: Stress over time. Lines are mean responses. Confidence intervals are within-confidence intervals.

We measure stress as the change in stress rating from pre to post induction. The following plot shows this:

the_data = make_the_data(my_data.Q.stressed)[, outcome := "stress"]
attr(the_data,"N")

Table 3.1: Number of particicipants per condition by condition and stage

Drug	state.condition	pre.induction	post.induction
oxycodone	control	62	62
oxycodone	stress	63	63
placebo	control	63	63
placebo	stress	63	63

the_data %>% plot_pre_post() + scale_col_stress

Figure 3.2: Change in stress ratings from pre to post state induction. Each line represents an individual. Stress is reflected in the slope of the lines.

We can further summarize this information by showing the distribution of the slopes in the previous figure:

p = plot_stresseff_raw(the_data)
raw_contrasts.stress = 
  p$data %>% 
  .[, .(participant,Sex,Drug, stress_effect)] %>% 
  .[, Outcome := "Stressed (VAS)"]
p

Figure 3.3: Distribution of change in stress ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

The next plot shows the individual data points, where each data point is calculated as a pre-post difference, together with mean and standard errors.

plot_stresseff_raw2(the_data)

Figure 3.4: Mean and CIs of change in stress ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

To estimate the effect of stress induction on stress ratings, we use a hierarchical ordered logistic model regression. The dependent variables are pre and post stress induction ratings and the independent variable the state induction. We model the responses as main effects of Stress and Sex plus the state-induction X sex interaction, where the interaction term measures if the effect of stress induction differed between males and females:

stress_formula = response.c ~ induced.state + Sex + stage:induced.state:Sex + (1|session) + (1|participant)

stress_formula = response.c ~ induced.state + Sex + stage:induced.state:Sex + (1|session) + (1|participant) 

We estimate a standard hierarchical linear regression with Gaussian errors.

ffit = lme4::lmer(response ~ induced.state + Sex + stage:induced.state:Sex + (1|session) + (1|participant) , data = the_data)

To check the model, we first investigate model predictions. The following figure shows that the model predictions, and even the best linear predictor predict impossible negative values.

error_sd = 
  VarCorr(ffit) %>% attr("sc")
p1 = 
  data.table(yhat = predict(ffit)) %>% 
  .[, lower := yhat - error_sd] %>% 
  .[, upper := yhat + error_sd] %>% 
  .[, state.condition := the_data$state.condition] %>% 
  .[order(yhat)] %>% 
  .[, j := 1:.N] %>% 
  ggplot(aes(x = j, y = yhat, ymin = lower, ymax = upper,
             color = state.condition, fill = state.condition)) + 
   geom_hline(yintercept = 0, color = "red") +
  geom_linerange(alpha = .5) +
  geom_point(size = .4) +
  xlab("") +
  ylab("fitted + error") + 
  scale_col_stress
p2 = 
  rbind(the_data[, .(response)][, data := "observed"],
           data.table(response = predict(ffit), data = "predicted")) %>% 
  ggplot(aes(x = response, fill = data)) + 
  geom_histogram(position = "identity", alpha = .5)
  
(p1 | p2) + plot_layout(widths = c(2,1))

Figure 3.5: Predictions of the linear model +/- 1sd of the normal error distribution.

Consistently, the following figure shows that the residuals are correlated with the outcome variable and not normally distributed.

res = data.table(residuals = residuals(ffit), fitted = predict(ffit))
p1 = res %>% 
  ggplot(aes(x = fitted, y = residuals)) + 
  geom_point() + 
  geom_smooth()
p2 = 
  res %>% 
  ggplot(aes(x = residuals)) + 
  geom_histogram(aes(y =..density..),
                 breaks = seq(-60, 60, by = 5), 
                 colour = "black", 
                 fill = "white") +
  stat_function(fun = dnorm,
                args = list(mean = mean(res$residuals),
                            sd = sd(res$residuals)))

p1 | p2

Figure 3.6: Residuals from a linear regression model. Left: Predicted values and residuals are not independent. Right: Residuals are not normally distributed.

rm(error_sd,ffit,p1,p2,res)

In sum, a linear regression model with Gaussian errors is not a good model for the rating scale data. We therefore use ordinal regression models for coarsened data. Specifically, we coarsen the questionnaire data into 20 equally sized bins and perform Bayesian ordinal logistic regressions on the coarsened data. To report results, we back-calculate effects (including credible intervals) of independent variables on the original 0-100 rating scale.

3.1.2 Ordinal logistic regression

We coarsen the outcome variable in preparation of the ordinal regression model

the_data = the_data %>% 
  coarsen(var = "response", bins = 20, range = c(0,100))

Here we run the Bayesian ordinal logistic regression for the first stress induction:

fn = "brmsP1/stressed_acat.Rdata"
if (file.exists(fn)) {
  load(fn)
} else {
  fit = fit_ordinal_model(the_data, stress_formula, fn, family = acat())
}

Click to show sample diagnostics

Table 3.2: Summary of model paramters

parameter	mean	sd	q5	q95	rhat	ess_bulk	ess_tail
Threshold[1]	1.3	1.13	-0.6	3.1	1.001	1729	2599
Threshold[2]	0.5	1.13	-1.4	2.3	1.001	1673	2187
Threshold[3]	0.4	1.15	-1.5	2.3	1.001	1765	2496
Threshold[4]	1.1	1.17	-0.9	3.0	1.002	1674	2509
Threshold[5]	0.3	1.18	-1.7	2.2	1.001	1823	2714
Threshold[6]	0.8	1.20	-1.2	2.7	1.001	1817	2730
Threshold[7]	0.6	1.23	-1.5	2.6	1.000	1809	2527
Threshold[8]	0.7	1.23	-1.3	2.7	1.001	1880	2510
Threshold[9]	0.4	1.22	-1.6	2.4	1.001	1831	2586
Threshold[10]	0.5	1.22	-1.4	2.5	1.001	1832	2726
Threshold[11]	0.8	1.22	-1.2	2.8	1.001	1797	2584
Threshold[12]	1.1	1.22	-1.0	3.1	1.001	1906	2855
Threshold[13]	0.8	1.25	-1.3	2.9	1.002	1903	2662
Threshold[14]	0.5	1.24	-1.6	2.6	1.001	1890	2659
Threshold[15]	0.8	1.22	-1.3	2.8	1.002	1826	2567
Threshold[16]	1.5	1.26	-0.5	3.6	1.001	1957	2567
Threshold[17]	1.2	1.29	-0.9	3.3	1.001	1998	2616
Threshold[18]	1.0	1.31	-1.2	3.2	1.001	1944	2541
Threshold[19]	0.7	1.30	-1.5	2.8	1.001	1975	2624
b_induced.statestress	0.4	1.33	-1.8	2.6	1.001	2033	2460
b_Sexwoman	0.1	1.23	-2.0	2.0	1.002	1956	2470
b_induced.statecontrol:Sexman:stagepre.induction	-0.2	1.11	-2.0	1.6	1.001	1695	2674
b_induced.statestress:Sexman:stagepre.induction	0.0	2.01	-3.4	3.3	1.000	6573	2797
b_induced.statecontrol:Sexwoman:stagepre.induction	0.0	1.14	-1.9	1.8	1.003	1612	2288
b_induced.statestress:Sexwoman:stagepre.induction	0.0	2.06	-3.4	3.4	1.000	5690	2653
b_induced.statecontrol:Sexman:stagepost.induction	-0.2	1.11	-2.0	1.6	1.001	1702	2693
b_induced.statestress:Sexman:stagepost.induction	0.1	1.34	-2.1	2.4	1.001	2249	2742
b_induced.statecontrol:Sexwoman:stagepost.induction	-0.1	1.14	-1.9	1.8	1.002	1611	2316
b_induced.statestress:Sexwoman:stagepost.induction	0.3	1.37	-1.9	2.5	1.003	2147	2807
sd_participant__Intercept	0.3	0.05	0.2	0.4	1.002	1429	2363
sd_session__Intercept	0.0	0.07	0.0	0.2	1.005	1290	1338
model: response.c ~ induced.state + Sex + stage:induced.state:Sex + (1 | session) + (1 | participant) b_ = fix effect parameters sd_ standard deviation of hierarchical parameters ess = effective sample size

The next figure shows that the model describes the observed data well.

#pp_check(fit, ndraws = 4000, type = "bars")

Due to non-collapsibility of OR we cannot simply interpret the interaction effect. Instead, We calculate the contrasts manually from posterior predictions:

tmp = plot_prepost_contrast(fit, the_data,"stress", my_contrast = "stress")
tmp[[1]]

Figure 3.7: Posterior distributions of estimated stress induction effects on self reported stress on the 0-100 VAS scale in women, men, and sex difference (women-men). Density plots show posterior distributions. Dots and lines identify means, and 50 & 95% credible intervals.

tmp$table

Table 3.3: Modeled post minus pre state induction on stress ratings by condition (control, stress) and the stress effect by sex.

Sex	control	stress	stress_effect
men	0 [-2, 3];	29 [24, 34];	29 [24, 34];
women	-0 [-3, 3];	47 [42, 52];	48 [42, 53];
total	-0 [-2, 2];	38 [34, 42];	38 [34, 42];
Sex difference	-1 [-4, 3];	18 [10, 25];	18 [11, 26];

post = tmp[[2]]
summary_table.stress = tmp$summary_table[, Outcome := "Stressed (VAS)"]

The estimated effect of stress from the state induction is NaN (CI=[NA, NA]; BFNA; P(S>C)NA). The effect is 48 (CI=[42, 53]; BF>10^5; P(S>C)>0.999) in females, 29 (CI=[24, 34]; BF>10^5; P(S>C)>0.999) in males and the sex difference is NaN (CI=[NA, NA]; BFNA; P(S>C)NA).

3.2 Negative mood

We repeat the same analysis as for stress and also take into account that there are multiple items measuring negative mood. These items are: distressed, anxious, vulnerable, embarassed, irritable

my_data.Q.neg_mood = select_data(negaffect_items)
my_data.Q.neg_mood %>% plot_stages()

Figure 3.8: Negative mood over time. Lines are mean responses.

the_data = make_the_data(my_data.Q.neg_mood)[, outcome := "neg_mood"]
attr(the_data,"N")

Table 3.4: Number of particicipants per condition by condition and stage

Drug	state.condition	pre.induction	post.induction
oxycodone	control	62	62
oxycodone	stress	63	63
placebo	control	63	63
placebo	stress	63	63

the_data %>% plot_pre_post() + scale_col_stress

Figure 3.9: Change in negative mood ratings from pre to post state induction. Each line represents an individual. Stress is reflected in the slope of the lines.

We can further summarize this information by showing the distribution of the slopes in the previous figure:

plot_stresseff_raw(the_data)

Figure 3.10: Distribution on change in negative mood ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

plot_stresseff_raw2(the_data)

Figure 3.11: Mean and CIs of change in negative mood ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

We extend the analysis model to account for multiple items.

stress_formula = response.c ~ stage:induced.state:Sex + session + (1|participant) + (1 | item_name)

stress_formula = response.c ~ stage:induced.state:Sex + (1|session) + (1|participant) + (1 | item_name)

Here we run the Bayesian ordinal logistic regression:

fn = "brmsP1/neg_mood_acat.Rdata"
the_data = the_data %>% 
  coarsen(var = "response", bins = 20, range = c(0,100))
if (file.exists(fn)) {
  load(fn)
} else {
  fit = fit_ordinal_model(the_data, stress_formula, fn, family = acat())
}
#pp_check(fit, ndraws = 2000, type = "bars")

Click to show sample diagnostics

Table 3.5: Summary of model paramters

parameter	mean	sd	q5	q95	rhat	ess_bulk	ess_tail
Threshold[1]	1.2	0.84	-0.2	2.5	1.003	1013	1426
Threshold[2]	0.8	0.84	-0.6	2.1	1.003	1032	1445
Threshold[3]	0.5	0.85	-0.9	1.9	1.003	1052	1352
Threshold[4]	0.6	0.86	-0.9	1.9	1.003	1045	1482
Threshold[5]	0.4	0.87	-1.0	1.8	1.003	1100	1510
Threshold[6]	0.2	0.88	-1.1	1.6	1.003	1084	1523
Threshold[7]	0.9	0.89	-0.6	2.3	1.003	1099	1567
Threshold[8]	0.5	0.92	-1.0	2.0	1.004	1134	1653
Threshold[9]	-0.3	0.90	-1.8	1.1	1.003	1150	1534
Threshold[10]	0.8	0.89	-0.7	2.2	1.003	1144	1593
Threshold[11]	1.0	0.92	-0.4	2.5	1.002	1151	1868
Threshold[12]	0.0	0.92	-1.5	1.5	1.003	1163	1557
Threshold[13]	0.4	0.90	-1.1	1.9	1.003	1150	1621
Threshold[14]	1.3	0.92	-0.2	2.8	1.003	1148	1795
Threshold[15]	0.7	0.95	-0.9	2.3	1.003	1215	1731
Threshold[16]	0.5	0.96	-1.1	2.1	1.004	1265	1964
Threshold[17]	2.4	1.11	0.6	4.2	1.001	1541	2185
Threshold[18]	1.5	1.33	-0.7	3.7	1.001	1815	2153
Threshold[19]	-1.8	1.25	-4.0	0.1	1.002	1672	2282
b_stagepre.induction:induced.statecontrol:Sexman	-0.3	0.82	-1.7	1.0	1.003	1060	1380
b_stagepost.induction:induced.statecontrol:Sexman	-0.4	0.82	-1.8	0.9	1.002	1075	1408
b_stagepre.induction:induced.statestress:Sexman	0.0	1.98	-3.3	3.2	1.001	6671	2940
b_stagepost.induction:induced.statestress:Sexman	0.2	0.82	-1.1	1.5	1.003	1060	1341
b_stagepre.induction:induced.statecontrol:Sexwoman	-0.1	0.82	-1.4	1.2	1.003	1062	1376
b_stagepost.induction:induced.statecontrol:Sexwoman	-0.2	0.82	-1.5	1.1	1.003	1060	1397
b_stagepre.induction:induced.statestress:Sexwoman	0.0	2.03	-3.3	3.4	1.001	6800	2443
b_stagepost.induction:induced.statestress:Sexwoman	0.4	0.82	-0.9	1.7	1.003	1062	1454
sd_item_name__Intercept	0.3	0.22	0.1	0.7	1.002	1180	1233
sd_participant__Intercept	0.3	0.04	0.2	0.3	1.002	1208	1921
sd_session__Intercept	0.1	0.11	0.0	0.3	1.005	1064	1478
model: response.c ~ stage:induced.state:Sex + (1 | session) + (1 | participant) + (1 | item_name) b_ = fix effect parameters sd_ standard deviation of hierarchical parameters ess = effective sample size

We calculate the contrasts manually from posterior predictions:

tmp = plot_prepost_contrast(fit, the_data,"neg_mood", my_contrast = "stress")
tmp[[1]]

Figure 3.12: Posterior distributions of estimated stress induction effects on self reported negative mood on the 0-100 VAS scale in women, men, and sex difference (women-men). Density plots show posterior distributions. Dots and lines identify means, and 50 & 95% credible intervals.

tmp$table

Table 3.6: Modeled post minus pre state induction on neg_mood ratings by condition (control, stress) and the stress effect by sex.

Sex	control	stress	stress_effect
men	-0 [-1, 0];	13 [11, 15];	13 [11, 16];
women	-1 [-2, 1];	25 [22, 28];	26 [23, 28];
total	-1 [-1, 0];	19 [17, 21];	20 [18, 21];
Sex difference	-0 [-2, 1];	12 [9, 15];	12 [9, 16];

post = tmp[[2]]
summary_table.stress = 
  rbind(summary_table.stress,
        tmp$summary_table[, Outcome := "Negative Mood"])

The estimated effect of stress from the state induction is NaN (CI=[NA, NA]; BFNA; P(S>C)NA). The effect is 26 (CI=[23, 28]; BF>10^5; P(S>C)>0.999) in females, 13 (CI=[11, 16]; BF>10^5; P(S>C)>0.999) in males and the sex difference is NaN (CI=[NA, NA]; BFNA; P(S>C)NA).

3.3 Positive mood

We repeat the same analysis as for stress and also take into account that there are multilple items measuring negative mood. These items are: distressed, anxious, vulnerable, embarassed, irritable

my_data.Q.pos_mood = select_data(posaffect_items)
my_data.Q.pos_mood %>% plot_stages()

Figure 3.13: Positive mood over time. Lines are mean responses.

the_data = make_the_data(my_data.Q.pos_mood)[, outcome := "pos_mood"]
attr(the_data,"N")

Table 3.7: Number of particicipants per condition by condition and stage

Drug	state.condition	pre.induction	post.induction
oxycodone	control	62	62
oxycodone	stress	63	63
placebo	control	63	63
placebo	stress	63	63

the_data %>% plot_pre_post() + scale_col_stress

Figure 3.14: Change in positive mood ratings from pre to post state induction. Each line represents an individual. The stress induction effect is reflected in the slope of the lines.

We can further summarize this information by showing the distribution of the slopes in the previous figure:

plot_stresseff_raw(the_data)

Figure 3.15: Distribution on change in positive mood ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

plot_stresseff_raw2(the_data)

Figure 3.16: Mean and CIs of change in positive mood ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Here we run the Bayesian ordinal logistic regression:

fn = "brmsP1/pos_mood_acat.Rdata"
the_data = the_data %>% 
  coarsen(var = "response", bins = 20, range = c(0,100))
if (file.exists(fn)) {
  load(fn)
} else {
  fit = fit_ordinal_model(the_data, stress_formula, fn, family = acat())
}
#pp_check(fit, ndraws = 2000, type = "bars")

Click to show sample diagnostics

Table 3.8: Summary of model paramters

parameter	mean	sd	q5	q95	rhat	ess_bulk	ess_tail
Threshold[1]	0.6	0.89	-0.8	2.1	1.010	659	1230
Threshold[2]	-0.5	0.87	-1.9	1.0	1.007	687	1079
Threshold[3]	-0.9	0.89	-2.3	0.6	1.009	681	1263
Threshold[4]	0.0	0.87	-1.4	1.5	1.009	662	1224
Threshold[5]	-0.5	0.87	-1.9	1.0	1.008	676	945
Threshold[6]	-0.5	0.88	-1.8	1.0	1.009	672	1219
Threshold[7]	-0.3	0.86	-1.6	1.2	1.010	645	950
Threshold[8]	-0.1	0.87	-1.5	1.3	1.008	670	1101
Threshold[9]	-0.9	0.86	-2.2	0.6	1.009	645	1094
Threshold[10]	0.0	0.86	-1.3	1.4	1.009	651	1201
Threshold[11]	0.0	0.86	-1.4	1.4	1.009	651	1155
Threshold[12]	-0.3	0.85	-1.6	1.2	1.009	645	1059
Threshold[13]	0.1	0.86	-1.3	1.5	1.009	643	1071
Threshold[14]	-0.1	0.86	-1.4	1.4	1.009	649	1105
Threshold[15]	0.1	0.85	-1.2	1.5	1.009	643	968
Threshold[16]	0.2	0.85	-1.1	1.7	1.010	637	1203
Threshold[17]	0.2	0.85	-1.2	1.6	1.009	647	1125
Threshold[18]	0.7	0.86	-0.6	2.1	1.009	654	1067
Threshold[19]	-0.7	0.86	-2.1	0.7	1.009	640	1150
b_stagepre.induction:induced.statecontrol:Sexman	0.2	0.84	-1.1	1.6	1.011	645	1059
b_stagepost.induction:induced.statecontrol:Sexman	0.2	0.84	-1.1	1.6	1.011	645	1021
b_stagepre.induction:induced.statestress:Sexman	0.0	2.03	-3.4	3.3	1.002	3811	2923
b_stagepost.induction:induced.statestress:Sexman	-0.1	0.84	-1.4	1.3	1.011	644	1088
b_stagepre.induction:induced.statecontrol:Sexwoman	0.1	0.84	-1.2	1.5	1.010	650	1136
b_stagepost.induction:induced.statecontrol:Sexwoman	0.1	0.84	-1.2	1.5	1.010	649	1105
b_stagepre.induction:induced.statestress:Sexwoman	-0.1	1.96	-3.3	3.2	1.000	3766	2717
b_stagepost.induction:induced.statestress:Sexwoman	-0.3	0.84	-1.6	1.1	1.010	650	1144
sd_item_name__Intercept	0.2	0.10	0.1	0.3	1.004	820	1189
sd_participant__Intercept	0.2	0.02	0.2	0.2	1.008	548	1200
sd_session__Intercept	0.0	0.06	0.0	0.1	1.001	768	1024
model: response.c ~ stage:induced.state:Sex + (1 | session) + (1 | participant) + (1 | item_name) b_ = fix effect parameters sd_ standard deviation of hierarchical parameters ess = effective sample size

We calculate the contrasts manually from posterior predictions:

tmp = plot_prepost_contrast(fit, the_data, "pos_mood", my_contrast = "stress")
tmp[[1]]

Figure 3.17: Posterior distributions of estimated stress induction effects on self reported positive mood on the 0-100 VAS scale in women, men, and sex difference (women-men). Density plots show posterior distributions. Dots and lines identify means, and 50 & 95% credible intervals.

tmp$table

Table 3.9: Modeled post minus pre state induction on pos_mood ratings by condition (control, stress) and the stress effect by sex.

Sex	control	stress	stress_effect
men	1 [-1, 3];	-19 [-21, -16];	-20 [-23, -17];
women	2 [-1, 4];	-33 [-36, -30];	-35 [-38, -32];
total	1 [-0, 3];	-26 [-28, -24];	-27 [-29, -25];
Sex difference	1 [-3, 4];	-14 [-18, -11];	-15 [-19, -11];

post = tmp[[2]]
summary_table.stress = 
  rbind(summary_table.stress,
        tmp$summary_table[, Outcome := "Positive Mood"])

The estimated effect of stress from the state induction is NaN (CI=[NA, NA]; BFNA; P(S>C)NA). The effect is -35 (CI=[-38, -32]; BF<10^-5; P(S>C)<0.001) in females, -20 (CI=[-23, -17]; BF<10^-5; P(S>C)<0.001) in males and the sex difference is NaN (CI=[NA, NA]; BFNA; P(S>C)NA).

3.4 Heart rate

See the script “HRV.R” for code for pre-processing of the raw heart rate data.

In the data.table HR the column time indicates seconds from the middle of the stress induction and BPS indicates heart beats per second. We down sample to measure heart rate only every 15s.

load("data/HRV/HRV.Rdata")
HR = 
  HR %>% 
  .[, participant := as.numeric(gsub("E","",Participant))] %>% 
  .[participant %in% unique(my_data$participant)] %>% 
  .[, state.condition := ifelse(Stress == "Stress", "stress", "control")] %>% 
  .[, Sex := ifelse(Gender == "Male","men","women")] %>% 
  .[, minute := (time/60)]  

tHR = 
  HR %>% 
  .[, time.c := cut(time,
                    breaks = seq(-1200,1200,by = 20),
                    include.lowest = TRUE)] %>% 
  .[, .(time = mean(time), BPS = mean(BPS)),
    by = .(Participant, Session,Sex, Drug, Condition, state.condition, time.c)] %>% 
  .[, time.c := NULL]

tmp = 
  HR %>% 
  .[abs(minute) < 3.5] %>% # stress induction period
  .[, .(mBPS = mean(BPS)), by = .(participant, Sex, Drug, state.condition)] %>% 
  dcast(participant + Sex + Drug ~ state.condition, value.var = "mBPS") %>% 
  .[, stress_effect := stress-control] %>% 
  .[, .(participant, Sex, Drug, stress_effect)] %>% 
  .[, Outcome := "Heart rate, BPM"]

raw_contrasts.stress = 
  rbind(
    raw_contrasts.stress,
    tmp
  )
rm(tmp)

To estimate the effect of stress we fit a Bayesian spline model and compare the estimated heart rate after removing session effects.

Before the model can be estimated, we prepare the date by

• re-scaling time and heart rate (BPS) to facilitate parameter estimation.
• creating a group variable grp by crossing the Stress, Drug, and Sex variable. This will allow us to fit splines for these “groups” and makes it possible to estimate stress effects stratified by drug condition and sex.
• creating a auto-regression group variable AR.grp by crossing the variables participant and session. This variable will be used to group observations for the estimation of autoregressive error terms.

m.BPS = mean(tHR$BPS)
sd.BPS = mean(tHR$BPS)
tHR %>% 
  .[, sBPS := (BPS-m.BPS)/sd.BPS] %>% 
  .[, time := time/1000] 

tHR %>% 
  .[, grp := paste(c(as.character(state.condition),Drug,Sex), collapse = ":"), by = 1:nrow(tHR)] %>% 
  .[, AR.grp := paste(c(Participant,Session), collapse = ":"), by = 1:nrow(tHR)]

The brms model is then specified as:

HRV.formula = sBPS ~ 
  s(time, by = grp) + # smooth effects stratified by sex, drug, and stress condition
  (1|Session) +
  (1|Participant) + # participant and session level random effects
  ar(time = time, gr = AR.grp, p = 1, cov = FALSE) # autoregressive error terms

This model includes an adjustment for session, nested in participants.

Next we estimate the model:

fn = "brmsP1/HRV_induction_2.Rdata"

if (file.exists(fn)) {
  load(fn)
} else {
  bfit =  
    brm(HRV.formula,
        family = gaussian(),
        backend = "cmdstanr",
        iter = 1000, 
        data = tHR,
        prior = prior(normal(0,2), class = "b"))
  save(bfit, file = fn)
}

Click to show sample diagnostics

Table 3.10: Summary of model paramters

parameter	mean	sd	q5	q95	rhat	ess_bulk	ess_tail
b_Intercept	0.0	0.02	0.0	0.0	1.027	168	366
bs_stime:grpcontrol:Oxycodone:man_1	0.0	0.40	-0.6	0.7	1.004	1211	1347
bs_stime:grpcontrol:Oxycodone:woman_1	-0.3	0.38	-0.9	0.3	1.003	1092	1259
bs_stime:grpcontrol:Placebo:man_1	0.0	0.41	-0.6	0.7	1.001	859	1209
bs_stime:grpcontrol:Placebo:woman_1	0.4	0.35	-0.2	0.9	1.002	1321	1267
bs_stime:grpstress:Oxycodone:man_1	-0.1	0.42	-0.8	0.6	1.004	1259	1202
bs_stime:grpstress:Oxycodone:woman_1	0.2	0.41	-0.5	0.8	1.000	1341	1391
bs_stime:grpstress:Placebo:man_1	0.8	0.38	0.2	1.5	1.001	1060	1395
bs_stime:grpstress:Placebo:woman_1	0.0	0.40	-0.6	0.7	1.001	1404	1268
sd_Participant__Intercept	0.1	0.01	0.1	0.1	1.029	227	411
sd_Session__Intercept	0.0	0.02	0.0	0.0	1.003	447	488
sigma	0.1	0.00	0.1	0.1	1.007	5462	1202
model: sBPS ~ s(time, by = grp) + (1 | Session) + (1 | Participant) + ar(time = time, gr = AR.grp, p = 1, cov = FALSE) b_ = fix effect parameters sd_ standard deviation of hierarchical parameters ess = effective sample size

Here are the expected heart rates split by state and drug conditions and sex:

fn = "brmsP1/HRV_induction_2_pp.Rdata"

# new data for model predictions, where all participants have the same times
new.data = unique(tHR[, .(grp,Participant, Session, AR.grp)])
times = seq(min(bfit$data$time), max(bfit$data$time), length.out = 51)
new.data = 
  do.call(rbind,
          lapply(1:nrow(new.data), 
                 function(x) cbind(new.data[x],time = times)))

# generate posterior predictions
if (file.exists(fn)) {
  load(fn)
} else {
  pp = posterior_epred(bfit, newdata = new.data, cores = 4)
  
  # merge predictions with  drug, stress, sex and time info
  pdata = cbind(new.data, t(pp)) %>% 
    melt(id.vars = names(new.data), variable.name = "iter") %>% 
    .[, iter := gsub("V","",iter)]
  
  # average over participants with the same drug, stress, sex values
  pdata = pdata %>% 
    .[, .(value = mean(value)), by = .(grp,iter, time)]
  
  # prepare data for plotting
  pdata %>% 
    .[, c("Stress","Drug","Gender") := tstrsplit(grp,":")] %>% 
    .[, state.condition := factor(Stress, levels = c("stress","control"))] %>% 
    .[, Drug := factor(Drug)] %>% 
    .[, Sex := factor(Gender)] %>% 
    .[, iter := factor(iter)] %>% 
    # bring outcome variable back to original scale
    .[, value := value*sd.BPS+m.BPS] 
  
  save(new.data, pp , pdata, file = fn)
}


# create a baseline corrected version of the outcome,
# where the drug/stress/sex "groups" have the same
# average heart rate in the minutes -20 to -10
pdata = 
  pdata %>% 
  .[, minute := (time*1000/60)] %>% 
  .[, bl := mean(value[minute < -10]), by = .(state.condition, Drug, Gender, iter)] %>% 
  .[, value.blc := value-bl+m.BPS] %>% 
  .[, sc.i := paste(state.condition,iter)]

pdata[, Sex := gsub("man","men",Sex)]

pdata = rbind(
  pdata,
  avg = pdata %>% 
    .[,
      .(value = fmean(value), bl = fmean(bl), value.blc = fmean(value.blc)),
      by = .(Stress,Drug,state.condition,time,minute,iter,sc.i)] %>% 
    .[, Sex := "avg"],
  fill = TRUE)

# plot the expected heart rates from the first 100 iterations
pdata %>% 
  .[iter %in% 1:100] %>% 
  .[, alpha := 1] %>% 
  ggplot(aes(x = minute, y = value.blc, color = state.condition, group = sc.i, alpha = alpha)) + 
  geom_vline(xintercept = c(-3.5,3.5)) +
  geom_line() + 
  scale_alpha(range = c(0,.1),guide="none") + 
  facet_grid(Sex~Drug) +
  xlab("Time in experiment") + 
  ylab("Beats per minute") + 
  scale_x_continuous(breaks = c(-20,0,20), labels = c("BL3","state induction","dose1")) + 
  scale_col_stress + 
  theme(legend.position = "top")

Figure 3.18: Estimated heart rate as a function of drug and state condition and sex. The black vertical lines indicate the start and the end of the stress induction period.

Next we calculate the stress effects by comparing heart rates during the induction period:

test.data = 
  pdata %>% 
  .[abs(minute)<3.5, # limit to induction period
    .(value.blc = mean(value.blc)), # calculate average of baseline corrected HR
    by = .(state.condition, Drug, Sex, iter)] %>% 
  # transform to long format
  melt(id.vars = c("state.condition", "Drug", "Sex", "iter"), variable.name = "outcome") %>%   # transform to wide format, so that one column is stress and the other control
  dcast(iter + Drug + Sex ~ state.condition, value.var = "value") %>% 
  # calculate stress effect
  .[, Stress_effect := stress - control]

sexdiff = 
  test.data %>% 
  .[Sex != "avg"] %>% 
  dcast(iter + Drug ~ Sex, value.var = "Stress_effect") %>% 
  .[, Stress_effect := women-men] %>% 
  .[, Sex := "women-men"]

test.data = rbind(
  test.data, sexdiff,
  fill = TRUE
)

test.data %>% 
  ggplot(aes(x = Stress_effect, fill = Sex)) + 
  stat_halfeye(alpha = .5) + 
  facet_grid(Drug~Sex) + 
  geom_vline(xintercept = 0) +
  scale_col_sex + 
  ylab("") + gg_no_y_axis

Figure 3.19: Posterior distributions of the effect of stress induction on heart rate stratified by drug condition and sex. Density plots show posterior distributions. Dots and lines identify means, and 50 & 95% credible intervals.

stress.effect = test.data[, .(Stress_effect = mean(Stress_effect)), by = .(Sex,iter)]
# update table with stress effects
old.nms = c("men", "women", "avg", "women-men")
new.nms = c("men", "women", "total", "Sex difference")
tmp = stress.effect[, .(get_mci(Stress_effect, get.P = ifelse(grepl("-",Sex),TRUE,FALSE))), by = .(Sex)] %>% 
  dcast(.~Sex) %>% 
  setnames(old.nms,new.nms) %>% 
  .[,new.nms,with = FALSE]
summary_table.stress = rbind(
  summary_table.stress,
  cbind(Outcome = "Heart rate, BPM",tmp)
)
rm(tmp)

On average, stress increases heart rate during the state induction by 12.4 (CI=[9.9, 14.9]; BF>10^5; P(S>C)>0.999) by beats per minutes. The effect in females was 15.1 (CI=[11.6, 18.5]; BF>10^5; P(S>C)>0.999) and 9.8 (CI=[6.4, 13.4]; BF>10^5; P(S>C)>0.999) in males the effect of sex on the stress effect is 5.2 (CI=[0.3, 10.5]; BF=47.8; P(S>C)=0.98).

Lastly, we formally test if the stress effects differs between drug conditions or between sexes.

contrast.data = 
  rbind(
  test.data %>% 
  dcast(iter + Drug ~ Sex, value.var = "Stress_effect") %>% 
  .[, estimate := women-men] %>% 
  .[, .(estimate = mean(estimate)), by = .(iter)] %>% 
  .[, contrast := "Stress:Sex (women-men)"],
  test.data %>% 
  dcast(iter + Sex ~ Drug, value.var = "Stress_effect") %>% 
  .[, estimate := Oxycodone-Placebo] %>% 
  .[, .(estimate = mean(estimate)), by = .(iter)] %>% 
  .[, contrast := "Stress:Drug (Oxycodone-Placebo)"])

contrast.data %>% 
  ggplot(aes(x = estimate)) + 
  geom_vline(xintercept = 0) +
  stat_halfeye(alpha = .5) + 
  facet_wrap(~contrast, scales = "free_x") +
  xlab("Estimated interaction effect") + 
  ylab("") + gg_no_y_axis

Figure 3.20: Interaction effects. Left: No clear differences in the effect of state induction between drug conditions. Right: State induction increases heart rate in females more than in males.

The stress induced heart rate increase was 5.2 (CI=[0.3, 10.5]; BF=47.8; P(S>C)=0.98) greater in females compared to males.

3.5 Cortisol levels

Compared to the no-stress control induction, social stress induction will increase cortisol levels (after the stress tasks).

 my_data.cortisol = 
  my_data.cortisol %>% 
  .[, concentr := mean(result), by = .(Sample.ID)] %>% 
  .[, cv.concentr := sd(result)/concentr, by = .(Sample.ID)] %>% 
  .[, .(participant, Drug, Stress, stage, plate, concentr, gender, session, time)] %>% 
  unique() %>% 
  .[, l.concentr := log(concentr)] %>% 
  .[, Sex := ifelse(gender == "f","women","men")] %>%
  .[, induced.state := ifelse(Stress == "stress" & stage > 3, "stress", "control")] %>% 
  .[, l.bl := mean(l.concentr[stage %in% 2:3]), by = .(participant, Drug, Stress)] %>% 
  .[, l.bl.m := mean(l.concentr[stage %in% 2:3], na.rm = TRUE)] %>% 
  .[, l.concentr.blc := l.concentr-l.bl+l.bl.m] %>% 
  .[, bl := mean(concentr[stage %in% 2:3]), by = .(participant, Drug, Stress)] %>% 
  .[, bl.m := mean(concentr[stage %in% 2:3], na.rm = TRUE)] %>% 
  .[, concentr.blc := concentr-bl+bl.m] %>% 
  setnames("Stress","state.condition") 

my_data.cortisol %>% 
  .[, .(N = length(unique(participant))), by = .(Drug,state.condition,stage)] %>%
  dcast(Drug+state.condition~stage, value.var = "N") %>% 
  kable(caption = "Number of participants with cortisol data by drug, state.condition and stage.") %>% kable_styling(full_width = FALSE) %>% 
  kable_classic()

Table 3.11: Number of participants with cortisol data by drug, state.condition and stage.

Drug	state.condition	2	3	4	6	8	10
oxycodone	control	56	56	57	58	58	58
oxycodone	stress	57	60	57	56	57	58
placebo	control	60	60	59	60	58	59
placebo	stress	60	60	56	56	58	57

my_data.cortisol %>% 
  .[, response := concentr] %>% 
  plot_stages(my_stages = 1:10, ymin.min = NA) + 
  ylab("Cortisol concentration (mg/dl)")

Figure 3.21: Cortisol levels over time. Lines are mean responses.

my_data.cortisol %>% 
  .[, response := concentr.blc] %>% 
  plot_stages(my_stages = 1:10, ymin.min = NA) + 
  ylab("Cortisol concentration (mg/dl)")

Figure 3.22: Baseline corrected cortisol levels over time. Lines are mean responses.

3.5.1 By drug

my_data.cortisol %>% 
  plot_stages(my_stages = 1:10, by.Drug = TRUE, fill = "Drug") + 
  theme_classic() + 
  ylab("Plasma cortisol concentration (mg/dl)") 

Figure 3.23: Cortisol levels over time by drug. Lines are mean responses.

From here on we use cortisol data collected until and including the second reminder for plotting. We do not use baseline corrected data for the statistical, because this leads to more missing data where participants do not have samples from the baseline period. In the statistical model baseline correction is achieved by adding a random intercept per participant.

my_data.cortisol[, response := concentr] 
the_data = 
  make_the_data(
    my_data.cortisol[stage < 10], n.post.stages = 10)[, outcome := "cortisol"][, session := factor(session)]
the_data %>% 
  plot_pre_post(by.Sex = FALSE, phase = "drug:induction") + 
  ylim(-1, 30) + 
  ylab("Difference plasma ortisol concentration (mg/dl)")

Figure 3.24: Change in cortisol levels from pre to post state induction. Each line represents an individual. Stress is reflected in the slope of the lines.

the_data %>% plot_pre_post(by.Sex = TRUE) +  ylim(-1, 30) + 
  ylab("Difference plasma cortisol concentration (mg/dl)") + scale_col_stress

Figure 3.25: Change in cortisol levels from pre to post state induction by sex. Each line represents an individual. Stress is reflected in the slope of the lines.

p = plot_stresseff_raw(the_data, by.Sex = TRUE)
p + ggtitle("stress_effect = cortisol_in_stress-corisol_in_control")

Figure 3.26: Distribution of change in plasma cortisol levels from pre to post state induction, stratified by drug condition (panels) and sex (colors).

attr(p,"table") %>% 
  kable(digits = 1,
        caption = "Effects of stress by drug and sex directly calculated from raw data.") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic()

Table 3.12: Effects of stress by drug and sex directly calculated from raw data.

Drug	Sex	N	stress_effect	sd	se	ci
oxycodone	men	24	2.3	4.0	0.8	1.7
oxycodone	women	31	2.0	3.4	0.6	1.3
placebo	men	26	3.0	3.9	0.8	1.6
placebo	women	31	2.1	3.4	0.6	1.3

tmp = 
  p$data %>% 
  .[, .(participant,Sex,Drug, stress_effect)] %>% 
  .[, Outcome := "Plasma cortisol, ng/ML"] %>% 
  .[, stress_effect := stress_effect*10] # get on the correct scale
raw_contrasts.stress = rbind(
  raw_contrasts.stress,
  tmp
)
rm(tmp)
p

Figure 3.27: Distribution of change in plasma cortisol levels from pre to post state induction, stratified by drug condition (panels) and sex (colors).

plot_stresseff_raw2(the_data, by.Sex = TRUE) 

Figure 3.28: Mean and CIs of change of plasma cortisol levels (mg/dl) from pre to post state induction, stratified by drug condition and sex (colors).

3.5.2 Statistical analysis: Five post induction measures

To precisely estimate the effect of stress induction on cortisol levels we estimate a hierarchical linear regression model on the cortisol levels (not baseline corrected, until including the 2nd reminder), using the model

cortisol_formula.drug = 
  bf(concentr ~ 0 + Sex:state.condition:Drug:stage + s(time, by = Sex) + (1|session) + (1 | participant) + (1 | plate) )
cortisol_formula.drug = 
  bf(concentr ~ 0 + Sex:state.condition:Drug:stage + s(time, by = Sex) + (1|session) + (state.condition:Drug:stage | participant) + (1 | plate) )

This model estimates log cortisol levels for the 16 “experimental cells” that are generated by crossing the factors state induction, drug and sex and pre/post induction. In addition, the model estimates effect of time of day and random effects for participant, stage, and session).

bfn = "brmsP1/cortisol_drug_plate+.Rdata"
if (file.exists(bfn)) {
  load(bfn) 
} else {
  my_prior = c(
    prior(normal(10,2), class = "b"),
    prior(normal(0,5), class = "b", coef = "stime:Sexman_1"),
    prior(normal(0,5), class = "b", coef = "stime:Sexwoman_1"),
    prior(normal(0,5), class = "sd", coef = "Intercept", group = "participant"),
    prior(normal(0,2), class = "sd"),
    prior(normal(0,5), class = "sds"),
    prior(normal(0,3), class = "sigma")
  )
  the_data %>% 
    .[, session := factor(session)] %>% 
    .[, participant := factor(participant)]
  bfit = 
    brm(cortisol_formula.drug,
        prior = my_prior,
        family = student(),
        data = the_data,
        iter = 2000,
        threads = threading(2, grainsize = ceiling(926/2)),
        control = list(adapt_delta = .8),
        backend = "cmdstanr")
  save(bfit, file = bfn)
}

Click to show sample diagnostics

Table 3.13: Summary of model paramters

parameter	mean	sd	q5	q95	rhat	ess_bulk	ess_tail
b_Sexman:state.conditioncontrol:Drugoxycodone:stagepre.induction	9.6	0.67	8.5	10.7	1.001	1416	2434
b_Sexwoman:state.conditioncontrol:Drugoxycodone:stagepre.induction	9.0	0.64	7.9	10.0	1.005	1187	1997
b_Sexman:state.conditionstress:Drugoxycodone:stagepre.induction	9.4	0.76	8.2	10.7	1.002	1428	2556
b_Sexwoman:state.conditionstress:Drugoxycodone:stagepre.induction	8.7	0.73	7.5	9.9	1.004	1416	2352
b_Sexman:state.conditioncontrol:Drugplacebo:stagepre.induction	10.4	0.81	9.1	11.8	1.000	1659	2562
b_Sexwoman:state.conditioncontrol:Drugplacebo:stagepre.induction	9.5	0.81	8.2	10.9	1.004	1518	2044
b_Sexman:state.conditionstress:Drugplacebo:stagepre.induction	10.4	0.81	9.1	11.7	1.000	1626	2337
b_Sexwoman:state.conditionstress:Drugplacebo:stagepre.induction	11.2	0.78	9.9	12.5	1.004	1441	2346
b_Sexman:state.conditioncontrol:Drugoxycodone:stagepost.induction	8.4	0.60	7.4	9.4	1.001	1342	2170
b_Sexwoman:state.conditioncontrol:Drugoxycodone:stagepost.induction	8.1	0.58	7.1	9.0	1.007	1107	1826
b_Sexman:state.conditionstress:Drugoxycodone:stagepost.induction	10.7	0.76	9.4	11.9	1.003	1552	2342
b_Sexwoman:state.conditionstress:Drugoxycodone:stagepost.induction	9.7	0.73	8.6	10.9	1.003	1389	2353
b_Sexman:state.conditioncontrol:Drugplacebo:stagepost.induction	9.6	0.67	8.5	10.7	1.001	1705	2617
b_Sexwoman:state.conditioncontrol:Drugplacebo:stagepost.induction	8.4	0.64	7.4	9.5	1.003	1564	2223
b_Sexman:state.conditionstress:Drugplacebo:stagepost.induction	12.2	0.90	10.7	13.7	1.000	1865	1959
b_Sexwoman:state.conditionstress:Drugplacebo:stagepost.induction	12.1	0.87	10.6	13.5	1.003	1636	2724
bs_stime:Sexman_1	-3.9	5.29	-12.6	5.0	1.000	2995	2331
bs_stime:Sexwoman_1	-4.0	4.63	-11.4	3.7	1.000	5263	2849
sd_participant__Intercept	4.4	0.48	3.7	5.3	1.001	1677	2213
sd_participant__state.conditioncontrol:Drugoxycodone:stagepre.induction	0.7	0.48	0.1	1.6	1.006	301	579
sd_participant__state.conditionstress:Drugoxycodone:stagepre.induction	2.1	0.49	1.2	2.9	1.003	759	843
sd_participant__state.conditioncontrol:Drugplacebo:stagepre.induction	3.5	0.50	2.7	4.4	1.001	904	1609
sd_participant__state.conditionstress:Drugplacebo:stagepre.induction	2.8	0.55	2.0	3.7	1.002	788	1826
sd_participant__state.conditioncontrol:Drugoxycodone:stagepost.induction	0.6	0.42	0.0	1.4	1.003	299	594
sd_participant__state.conditionstress:Drugoxycodone:stagepost.induction	2.8	0.40	2.2	3.5	1.001	1045	1535
sd_participant__state.conditioncontrol:Drugplacebo:stagepost.induction	3.0	0.38	2.4	3.6	1.001	549	740
sd_participant__state.conditionstress:Drugplacebo:stagepost.induction	4.1	0.44	3.4	4.8	1.004	1094	1297
sd_plate__Intercept	1.3	0.22	0.9	1.6	1.001	2110	2790
sd_session__Intercept	0.5	0.42	0.0	1.3	1.000	1726	2263
sigma	1.0	0.09	0.9	1.2	1.003	728	1947
model: concentr ~ 0 + Sex:state.condition:Drug:stage + s(time, by = Sex) + (1 | session) + (state.condition:Drug:stage | participant) + (1 | plate) b_ = fix effect parameters sd_ standard deviation of hierarchical parameters ess = effective sample size

pdata = 
  conditional_effects(bfit,"time:Sex")$time %>% 
  data.table() %>% 
  .[, .(Sex,time,estimate__,lower__,upper__)] %>% 
  .[, h := time*24] %>% 
  .[, Sex := gsub("man","men",Sex)]
pdata %>% 
  ggplot(aes(x = h, y = estimate__, fill = Sex, color = Sex)) + 
  geom_ribbon(aes(ymin = lower__, ymax = upper__), alpha = .5, color = NA) + 
  geom_line(size = 1) + 
  ylab("Plasma cortisol level") + 
  xlab("Time") + 
  scale_col_sex +
  scale_x_continuous(breaks = seq(10,20,2),
                     labels = paste0(seq(10,20,2),":00"))

Figure 3.29: estimated effect of time

The following figure shows a posterior predictive check, which confirms that the model captures drug and stress effects well.

draws = 
  bfit$fit %>% 
  as_draws() %>% 
  subset_draws("b_Sex", regex = TRUE) %>% 
  as_draws_matrix() %>% 
  data.table() %>% 
  .[, iter := 1:.N] %>% 
  melt(id.vars = "iter", variable.name = "cond", value.name = "est") %>% 
  .[grepl("pre|post",cond)] %>% 
  .[, cond := gsub("b_Sex|Drug|state.condition|stage","",cond)] %>% 
  .[, c("Sex","state.condition","Drug","stage") := tstrsplit(cond,":")] %>% 
  .[, cond := NULL] %>% 
  .[, Sex := gsub("man","men",Sex)]


obs = the_data[, .(est = mean(concentr)), by = c("Sex","state.condition","Drug","stage")]


obs %>% 
  ggplot(aes(x = stage, y=est, color = state.condition)) + 
  geom_point(aes(shape = state.condition),position = position_dodge(.2)) +
  geom_line(aes(lty = state.condition,group = factor(state.condition):factor(Drug)),
            position = position_dodge(.2)) +
  facet_grid(Drug~Sex) + 
  stat_interval(position = position_dodge(.2), data = draws, alpha = .25) +
  facet_grid(Drug~Sex) +
  scale_col_stress + 
  ylab("Cortisol level")

Figure 3.30: Observed and estimated cortisol levels pre and post state induction. Points are means of observed data and confidence bands show the 95% credible intervals of the estimated cortisol levels, after adjustment for session and time effects.

We use the estimated cortisol levels in the 16 cells (2 drug conditions * 2 state induction conditions * 2 sexes * 2 times) to calculate the contrasts of interest.

1. We calculate the effect of induction as the post-pre differences in the eight cells
2. We calculate the stress effect by subtracting the post-pre difference in the control condition from the post-pre in the stress condition, separately for Drug condition and Sex
3. From these 4 stress effects, we can calculate average stress effects and comparisons between drug conditions or sexes

draws = 
  bfit$fit %>% 
  as_draws() %>% 
  subset_draws("b_Sex", regex = TRUE) %>% 
  as_draws_matrix() %>% 
  data.table()

setnames(draws,
         names(draws),
         gsub("b_Sex|Drug|state.condition|stage|.induction","",names(draws)))
setnames(draws,
         names(draws),
         gsub(":",".",names(draws)))
setnames(draws,names(draws),gsub("man","men", names(draws)))

make_contrasts = function(draws) {
  # put onto scale used for paper
  draws = draws * 10
  tmp = draws %>% 
    .[, `:=`(
      men.stress.oxycodone = men.stress.oxycodone.post-men.stress.oxycodone.pre,
      men.stress.placebo = men.stress.placebo.post-men.stress.placebo.pre,
      men.control.oxycodone = men.control.oxycodone.post-men.control.oxycodone.pre,
      men.control.placebo = men.control.placebo.post-men.control.placebo.pre,
      women.stress.oxycodone = women.stress.oxycodone.post-women.stress.oxycodone.pre,
      women.stress.placebo = women.stress.placebo.post-women.stress.placebo.pre,
      women.control.oxycodone = women.control.oxycodone.post-women.control.oxycodone.pre,
      women.control.placebo = women.control.placebo.post-women.control.placebo.pre
    )] %>% 
    .[, .(
      men.stress_effect.oxycodone = men.stress.oxycodone - men.control.oxycodone,
      men.stress_effect.placebo = men.stress.placebo - men.control.placebo,
      women.stress_effect.oxycodone = women.stress.oxycodone - women.control.oxycodone,
      women.stress_effect.placebo = women.stress.placebo - women.control.placebo,
      men.drug_effect.stress = men.stress.oxycodone - men.stress.placebo,
      men.drug_effect.control = men.control.oxycodone - men.control.placebo,
      women.drug_effect.stress = women.stress.oxycodone - women.stress.placebo,
      women.drug_effect.control = women.control.oxycodone - women.control.placebo
    )] %>% 
    .[,`:=`(men.drug_x_stress = men.stress_effect.oxycodone - men.stress_effect.placebo,
            women.drug_x_stress = women.stress_effect.oxycodone - women.stress_effect.placebo
    )] %>% 
    .[, iter := 1:.N] %>% 
    melt(id.vars = "iter", value.name = "est", variable.name = "contr") %>% 
    .[, Sex := ifelse(grepl("women",contr),"women","men")] %>% 
    .[, contr := gsub("women.|men.","",contr)] %>% 
    dcast(iter + contr ~ Sex, value.var = "est") %>% 
    .[, `:=`(avg = (men+women)/2, `women-men` = women-men)] %>% 
    melt(id.var = c("iter","contr"), value.name = "est", variable.name = "Sex") %>%
    dcast(iter + Sex ~ contr) %>% 
    .[, stress_effect := (stress_effect.oxycodone+stress_effect.placebo)/2] %>% 
    .[, drug_effect := (drug_effect.stress+drug_effect.control)/2]
  return(tmp)
}

contrast.samples = make_contrasts(draws)

Here are main and interaction effects of state induction and drug:

contrast.samples %>% 
  .[, .(iter,Sex,drug_effect,stress_effect,drug_x_stress)] %>% 
  melt(id.var = c("iter","Sex"), variable.name = "contrast", value.name = "est") %>% 
  .[Sex == "avg"] %>% 
  ggplot(aes(x = est)) +
  facet_wrap(~contrast, nrow = 1) +
  geom_vline(xintercept = 0) +
  stat_halfeye() + 
  coord_flip() + 
  xlab("change of cortisol level (ng/ML)") +
  ylab("Posterior density")

Figure 3.31: Stress and drugs effects on log cortisol levels

write.cortisol.stat = function(var) {
  contrast.samples[Sex == "avg",get(var)] %>% 
    print_effect(digits = 1, BF.thresh = 0)
}

The estimated effect of stress from the state induction on log cortisol levels in ng/ML is over all 22.9 (CI=[14.4, 31.7]; BF>10^5; P(S>C)>0.999) (cortisol is higher under stress), the effect of oxycodone is -1.7 (CI=[-9.5, 5.7]; BF=0.5; P(S>C)=0.32) (cortisol is not credibly lower under oxycodone), and the interaction effect is -1.1 (CI=[-15.3, 13.2]; BF=0.8; P(S>C)=0.44).

The next figure shows that there are no large Sex difference in the cortisol effects:

contrast.samples %>% 
  .[, .(iter,Sex,drug_effect,stress_effect,drug_x_stress)] %>% 
  melt(id.var = c("iter","Sex"), variable.name = "contrast", value.name = "est") %>% 
  .[Sex %in% c("men","women")] %>% 
  ggplot(aes(x = est, fill = Sex)) +
  facet_wrap(~contrast, nrow = 1) +
  geom_vline(xintercept = 0) +
  stat_halfeye(alpha = .5) + 
  coord_flip() + 
  xlab("change of cortisol level (ng/ML)") + 
  scale_col_sex +
  ylab("Posterior density")

Figure 3.32: Stress and drugs effects on log cortisol levels, stratified by sex.

Here is the table with the relevant statistics.

contrast.samples %>% 
  .[Sex == "women-men",
    .(iter,Sex,drug_effect,stress_effect,drug_x_stress)] %>% 
  melt(id.var = c("iter","Sex"), variable.name = "contrast", value.name = "sex.diff") %>% 
  .[, .(estimate = mean(sex.diff),
        CI_low = quantile(sex.diff, probs = .025),
        CI_high = quantile(sex.diff, probs = .975),
        `P > 0` = mean(sex.diff > 0)),
    by = contrast] %>% 
  kable(digits = 2,
        caption = "Sex difference in cortisol analyses. Numbers are sex difference for the contrasts in ng/ML.") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic_2() 

Table 3.14: Sex difference in cortisol analyses. Numbers are sex difference for the contrasts in ng/ML.

contrast	estimate	CI_low	CI_high	P > 0
drug_effect	6.15	-8.76	20.74	0.79
stress_effect	-6.84	-24.21	10.27	0.22
drug_x_stress	1.80	-26.69	29.97	0.55

# update table with stress effects
tmp = 
  contrast.samples[, .(Sex,stress_effect)] %>% 
  .[, .(get_mci(stress_effect, get.P = ifelse(grepl("-",Sex),TRUE,FALSE))), by = .(Sex)] %>% 
  dcast(.~Sex) %>% 
  setnames(old.nms,new.nms) %>% 
  .[,new.nms,with = FALSE]
summary_table.stress = rbind(
  summary_table.stress,
  cbind(Outcome = "Plasma cortisol, ng/ML",tmp)
)
rm(tmp)

3.6 Summary of stress effects

summary_table.stress %>% 
  kable(caption = "Effects of Stress") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic() %>% 
  add_footnote("Numbers are mean [Credible interval] and posterior probability that an effect is larger than 0. The latter is only given for sex differences.")

Table 3.15: Effects of Stress

Outcome	men	women	total	Sex difference
Stressed (VAS)	29 [24, 34];	48 [42, 53];	38 [34, 42];	18 [11, 26];
Negative Mood	13 [11, 16];	26 [23, 28];	20 [18, 21];	12 [9, 16];
Positive Mood	-20 [-23, -17];	-35 [-38, -32];	-27 [-29, -25];	-15 [-19, -11];
Heart rate, BPM	10 [6, 13];	15 [12, 19];	12 [10, 15];	5 [0, 11]; 0.9795
Plasma cortisol, ng/ML	26 [14, 39];	19 [7, 31];	23 [14, 32];	-7 [-24, 10]; 0.21875
a Numbers are mean [Credible interval] and posterior probability that an effect is larger than 0. The latter is only given for sex differences.

summary_stats.stress = 
  summary_table.stress %>% 
  .[, .(Outcome, men, women)] %>% 
  melt(id.var = "Outcome", variable.name = "Sex") %>% 
  .[, value := gsub("\\]|;","",value)] %>% 
  .[, c("stress_effect","lower","upper") := tstrsplit(value,", | \\[")] %>% 
  .[, value := NULL] %>% 
  .[, `:=`(stress_effect = as.numeric(stress_effect), lower = as.numeric(lower), upper = as.numeric(upper))] %>% 
  .[Outcome %in% unique(raw_contrasts.stress$Outcome)]
raw_contrasts.stress %>% 
  .[, Sex := gsub("a","e",Sex)] %>% 
  .[, l := fmean(stress_effect)-3*fsd(stress_effect), by = .(Outcome,Sex)] %>% 
  .[, u := fmean(stress_effect)+3*fsd(stress_effect), by = .(Outcome,Sex)] %>% 
  .[, plot := stress_effect < u & stress_effect > l] %>% 
  .[, Drug := stringr::str_to_title(Drug)]

summary_stats.stress[, Outcome := paste("\u0394 ", Outcome)]
raw_contrasts.stress[, Outcome := paste("\u0394 ", Outcome)]
breaks_stress_fun <- function(x) {
  s.count <<- s.count + 1L
  switch(
    s.count,
    seq(0,60,20), c(0),
    seq(-50,150,50),c(0),
    seq(0,100,25), c(0),
    seq(0,60,20), c(0),
    seq(-50,150,50),c(0),
    seq(0,100,25)
  )
}
s.count = 0
summary_plot.stress = 
  summary_stats.stress %>% 
  ggplot(aes(x = Sex, y = stress_effect)) + 
  geom_bar(stat = "identity", aes(fill = Sex)) + 
  geom_hline(yintercept = 0) +
  geom_errorbar(aes(ymin = lower, ymax = upper), 
                linewidth = 1, width = 1/5) +
  facet_grid(Outcome~., scale = "free_y",switch = "y") + 
  scale_y_continuous(breaks = breaks_stress_fun) +
  scale_shape_manual(values = c(3,4)) +
  geom_quasirandom(data = raw_contrasts.stress[stress_effect < u & stress_effect > l],
                   aes(shape = Drug), color = adjustcolor("black",alpha = .5)) +
  guides(fill = guide_legend(nrow = 1, title.position = "left"), 
#         shape = guide_legend(nrow = 1, title.position = "left"),
         shape = "none") + 
  theme(legend.position = "top", axis.line.y = element_line(), 
        panel.grid.major.y = element_blank(),panel.spacing = unit(1, "lines"),
        strip.background = element_blank(), legend.key.size = unit(1,"line"),
        strip.placement = "outside", legend.title=element_text(size=12)) + 
  scale_col_sex2 + ylab("") + xlab("") + gg_no_x_axis

3.7 Reminders

3.7.1 Embarrassed (1st reminder)

We repeat the same analysis as for stress for the first reminder.

my_data.Q.emb = select_data(c("embarassed"))
my_stages = 5:6
my_data.Q.emb %>% plot_stages(my_stages = my_stages)

Figure 3.33: Embarrassed ratings over time. Lines are mean responses.

the_data = make_the_data(my_data.Q.emb, phase = "reminder1")[, outcome := "emb"]
attr(the_data,"N")

Table 3.16: Number of particicipants per condition by condition and stage

Drug	state.condition	pre.reminder1	post.reminder1
oxycodone	control	62	62
oxycodone	stress	63	63
placebo	control	63	63
placebo	stress	63	63

the_data %>% plot_pre_post("reminder") + scale_col_stress

Figure 3.34: Change in embarrassed ratings from pre to post state induction. Each line represents an individual. The stress induction effect is reflected in the slope of the lines.

We can further summarize this information by showing the distribution of the slopes in the previous figure:

plot_stresseff_raw(the_data)

Figure 3.35: Distribution on change in embarrassed ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

plot_stresseff_raw2(the_data)

Figure 3.36: Mean and CIs of change in embarrassed ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Here we run the Bayesian ordinal logistic regression:

fn = "brmsP1/Reminder/emb1_acat.Rdata"
stress_formula = response.c ~ stage:induced.state:Sex + session + (1|participant) 
the_data = the_data %>% 
  coarsen(var = "response", bins = 20, range = c(0,100))
if (file.exists(fn)) {
  load(fn)
} else {
  fit = fit_ordinal_model(the_data, stress_formula, fn, family = acat())
}
#pp_check(fit, ndraws = 2000, type = "bars")

Click to show sample diagnostics

Table 3.17: Summary of model paramters

parameter	mean	sd	q5	q95	rhat	ess_bulk	ess_tail
Threshold[1]	1.6	0.71	0.4	2.8	1.003	1066	1637
Threshold[2]	0.7	0.72	-0.5	1.9	1.005	1102	1720
Threshold[3]	0.4	0.74	-0.8	1.6	1.002	1169	1770
Threshold[4]	0.7	0.77	-0.6	2.0	1.002	1198	1912
Threshold[5]	0.8	0.79	-0.5	2.0	1.005	1217	1977
Threshold[6]	0.9	0.82	-0.5	2.3	1.001	1317	2212
Threshold[7]	1.0	0.90	-0.4	2.6	1.002	1518	2177
Threshold[8]	0.8	0.95	-0.7	2.4	1.003	1556	2354
Threshold[9]	-0.3	0.91	-1.8	1.1	1.001	1601	2542
Threshold[10]	2.0	0.93	0.5	3.5	1.000	1615	2377
Threshold[11]	0.7	1.05	-1.0	2.4	1.001	1636	2216
Threshold[12]	0.0	0.99	-1.7	1.6	1.001	1951	2550
Threshold[13]	0.8	0.90	-0.7	2.3	1.002	1530	2381
Threshold[14]	0.5	0.87	-0.9	2.0	1.004	1442	2194
Threshold[15]	3.1	1.14	1.3	5.1	1.001	2045	2822
Threshold[16]	0.2	1.25	-2.0	2.2	1.001	2040	2630
Threshold[17]	0.5	1.08	-1.3	2.2	1.002	2027	2547
Threshold[18]	0.4	1.00	-1.3	2.0	1.000	1811	2436
b_session2	0.1	0.07	0.0	0.2	1.001	2641	3169
b_session3	0.1	0.05	0.0	0.2	1.000	4649	3033
b_session4	0.1	0.07	0.0	0.2	1.000	3124	3268
b_stagepre.reminder1:induced.statecontrol:Sexman	-0.2	0.70	-1.4	0.9	1.004	1034	1569
b_stagepost.reminder1:induced.statecontrol:Sexman	-0.4	0.70	-1.5	0.8	1.003	1038	1613
b_stagepre.reminder1:induced.statestress:Sexman	0.0	0.69	-1.2	1.1	1.004	1020	1515
b_stagepost.reminder1:induced.statestress:Sexman	0.2	0.69	-0.9	1.4	1.004	1025	1660
b_stagepre.reminder1:induced.statecontrol:Sexwoman	-0.1	0.70	-1.2	1.1	1.004	1028	1736
b_stagepost.reminder1:induced.statecontrol:Sexwoman	-0.5	0.70	-1.6	0.7	1.004	1029	1527
b_stagepre.reminder1:induced.statestress:Sexwoman	0.3	0.69	-0.8	1.4	1.005	1017	1652
b_stagepost.reminder1:induced.statestress:Sexwoman	0.6	0.69	-0.5	1.7	1.005	1016	1639
sd_participant__Intercept	0.3	0.06	0.2	0.5	1.002	1196	2300
model: response.c ~ stage:induced.state:Sex + session + (1 | participant) b_ = fix effect parameters sd_ standard deviation of hierarchical parameters ess = effective sample size

We calculate the contrasts manually from posterior predictions:

tmp = plot_prepost_contrast(fit, the_data,"emb1", my_contrast = "stress")
tmp[[1]]

Figure 3.37: Posterior distributions of estimated first stress reminder effects on self reported embarrassment on the 0-100 VAS scale in women, men, and sex difference (women-men). Density plots show posterior distributions. Dots and lines identify means, and 50 & 95% credible intervals.

tmp$table

Table 3.18: Modeled post minus pre state induction on emb1 ratings by condition (control, stress) and the stress effect by sex.

Sex	control	stress	stress_effect
men	-1 [-4, 1];	8 [3, 12];	9 [4, 14];
women	-3 [-6, -1];	20 [13, 26];	23 [16, 29];
total	-2 [-4, -0];	14 [10, 18];	16 [11, 20];
Sex difference	-2 [-5, 2];	12 [4, 19];	14 [5, 22];

post = tmp[[2]]

The estimated effect of stress from the first state reminder is NaN (CI=[NA, NA]; BFNA; P(S>C)NA). The effect is 23 (CI=[16, 29]; BF>10^5; P(S>C)>0.999) in females, 9 (CI=[4, 14]; BF=1999; P(S>C)>0.999) in males and the sex difference is NaN (CI=[NA, NA]; BFNA; P(S>C)NA).

3.7.2 Embarrassed

We repeat the same analysis as for stress and also take into account that there are multiple reminders.

my_data.Q.emb = select_data(c("embarassed"))
my_stages = 5:8
my_data.Q.emb %>% plot_stages(my_stages = my_stages)

Figure 3.38: Embarrassed ratings over time. Lines are mean responses.

the_data = make_the_data(my_data.Q.emb, phase = "reminder")[, outcome := "emb"]
attr(the_data,"N")

Table 3.19: Number of particicipants per condition by condition and stage

Drug	state.condition	pre.reminder	post.reminder
oxycodone	control	62	62
oxycodone	stress	63	63
placebo	control	63	63
placebo	stress	63	63

the_data %>% plot_pre_post("reminder") + scale_col_stress

Figure 3.39: Change in embarrassed ratings from pre to post state induction. Each line represents an individual. The stress induction effect is reflected in the slope of the lines.

We can further summarize this information by showing the distribution of the slopes in the previous figure:

plot_stresseff_raw(the_data)

Figure 3.40: Distribution on change in embarrassed ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

plot_stresseff_raw2(the_data)

Figure 3.41: Mean and CIs of change in embarrassed ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Here we run the Bayesian ordinal logistic regression:

fn = "brmsP1/Reminder/emb_acat.Rdata"
stress_formula = response.c ~ induced.state + Sex + induced.state:Sex + session + (1|participant) + (1|stage)
the_data = the_data %>% 
  coarsen(var = "response", bins = 20, range = c(0,100))
if (file.exists(fn)) {
  load(fn)
} else {
  fit = fit_ordinal_model(the_data, stress_formula, fn, family = acat())
}
#pp_check(fit, ndraws = 2000, type = "bars")

Click to show sample diagnostics

Table 3.20: Summary of model paramters

parameter	mean	sd	q5	q95	rhat	ess_bulk	ess_tail
Threshold[1]	1.8	0.28	1.2	2.2	1.015	427	349
Threshold[2]	0.8	0.30	0.3	1.3	1.013	458	340
Threshold[3]	0.7	0.33	0.2	1.2	1.009	519	400
Threshold[4]	0.5	0.34	-0.1	1.1	1.008	555	511
Threshold[5]	0.8	0.37	0.2	1.4	1.005	754	778
Threshold[6]	1.1	0.41	0.4	1.8	1.007	665	469
Threshold[7]	1.0	0.50	0.2	1.9	1.001	1109	1415
Threshold[8]	0.5	0.52	-0.4	1.3	1.004	1148	1678
Threshold[9]	0.4	0.47	-0.4	1.1	1.004	987	1474
Threshold[10]	0.9	0.48	0.1	1.6	1.003	789	497
Threshold[11]	0.5	0.47	-0.3	1.3	1.003	954	1305
Threshold[12]	0.7	0.46	0.0	1.5	1.003	1112	1955
Threshold[13]	1.4	0.54	0.5	2.3	1.003	1357	1638
Threshold[14]	0.3	0.55	-0.6	1.2	1.002	1334	1980
Threshold[15]	1.9	0.63	0.9	3.0	1.003	1459	1893
Threshold[16]	0.8	0.73	-0.4	2.0	1.003	884	531
Threshold[17]	0.6	0.71	-0.6	1.8	1.001	1840	1558
Threshold[18]	1.4	0.75	0.2	2.7	1.004	1152	975
Threshold[19]	-0.6	0.68	-1.8	0.4	1.000	1708	1889
b_induced.statestress	0.4	0.06	0.3	0.5	1.003	1550	951
b_Sexwoman	0.2	0.11	0.0	0.3	1.012	696	1291
b_session2	0.0	0.04	-0.1	0.1	1.003	1640	2245
b_session3	0.0	0.03	0.0	0.1	1.003	3359	2516
b_session4	-0.1	0.04	-0.2	0.0	1.002	1559	2418
b_induced.statestress:Sexwoman	0.1	0.07	0.0	0.2	1.002	2012	1613
sd_participant__Intercept	0.3	0.05	0.2	0.4	1.008	675	1675
sd_stage__Intercept	0.4	0.35	0.1	1.1	1.017	386	1724
model: response.c ~ induced.state + Sex + induced.state:Sex + session + (1 | participant) + (1 | stage) b_ = fix effect parameters sd_ standard deviation of hierarchical parameters ess = effective sample size

We calculate the contrasts manually from posterior predictions:

tmp = plot_prepost_contrast(fit, the_data,"emb", my_contrast = "stress")
tmp[[1]]

Figure 3.42: Posterior distributions of estimated stress reminders effects on self reported embarrassment on the 0-100 VAS scale in women, men, and sex difference (women-men). Density plots show posterior distributions. Dots and lines identify means, and 50 & 95% credible intervals.

tmp$table

Table 3.21: Modeled post minus pre state induction on emb ratings by condition (control, stress) and the stress effect by sex.

Sex	control	stress	stress_effect
men	1 [1, 2];	7 [5, 9];	5 [4, 8];
women	2 [1, 3];	14 [10, 18];	12 [9, 16];
total	2 [1, 2];	11 [8, 14];	9 [6, 12];
Sex difference	0 [-0, 1];	7 [5, 10];	7 [4, 10];

post = tmp[[2]]

The estimated effect of stress from the state reminder on feeling embarrased is NaN (CI=[NA, NA]; BFNA; P(S>C)NA). The effect is 12 (CI=[9, 16]; BF>10^5; P(S>C)>0.999) in females, 5 (CI=[4, 8]; BF>10^5; P(S>C)>0.999) in males and the sex difference is NaN (CI=[NA, NA]; BFNA; P(S>C)NA).

3.7.3 Selfconscious

We repeat the same analysis as for stress and also take into account that there are multiple reminders.

my_data.Q.selfc = select_data(c("selfconscious"))
my_stages = 5:8
my_data.Q.selfc %>% plot_stages(my_stages=my_stages)

Figure 3.43: Embarrassed ratings over time. Lines are mean responses.

the_data = make_the_data(my_data.Q.selfc, phase = "reminder")[, outcome := "selfc"]
attr(the_data,"N")

Table 3.22: Number of particicipants per condition by condition and stage

Drug	state.condition	pre.reminder	post.reminder
oxycodone	control	62	62
oxycodone	stress	63	63
placebo	control	63	63
placebo	stress	63	63

the_data %>% plot_pre_post("reminder") + scale_col_stress

Figure 3.44: Change in selfconscious ratings from pre to post state induction. Each line represents an individual. The stress induction effect is reflected in the slope of the lines.

We can further summarize this information by showing the distribution of the slopes in the previous figure:

plot_stresseff_raw(the_data)

Figure 3.45: Distribution on change in embarrassed ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

plot_stresseff_raw2(the_data)

Figure 3.46: Mean and CIs of change in selfconcious ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Here we run the Bayesian ordinal logistic regression:

fn = "brmsP1/Reminder/selfc_acat.Rdata"
stress_formula = response.c ~ induced.state + Sex + induced.state:Sex + session + (1|participant) + (1|stage)
the_data = the_data %>% 
  coarsen(var = "response", bins = 20, range = c(0,100))
if (file.exists(fn)) {
  load(fn)
} else {
  fit = fit_ordinal_model(the_data, stress_formula, fn, family = acat())
}
#pp_check(fit, ndraws = 2000, type = "bars")

We calculate the contrasts manually from posterior predictions:

tmp = plot_prepost_contrast(fit, the_data,"selfc", my_contrast = "stress")
tmp[[1]]

Figure 3.47: Stress effects in women, men, and sex difference (women-men)

tmp$table

Table 3.23: Modeled post minus pre state induction on selfc ratings by condition (control, stress) and the stress effect by sex.

Sex	control	stress	stress_effect
men	1 [-0, 3];	1 [-0, 3];	-0 [-0, 0];
women	2 [-0, 4];	2 [-0, 4];	-0 [-0, 0];
total	2 [-0, 4];	2 [-0, 4];	-0 [-0, 0];
Sex difference	0 [-0, 1];	0 [-0, 1];	0 [-0, 0];

post = tmp[[2]]

The estimated effect of stress from the state reminder on feeling selfconscious is NaN (CI=[NA, NA]; BFNA; P(S>C)NA). The effect is 0 (CI=[0, 0]; BF=0.1; P(S>C)=0.06) in females, 0 (CI=[0, 0]; BF=0.1; P(S>C)=0.06) in males and the sex difference is NaN (CI=[NA, NA]; BFNA; P(S>C)NA).

3.7.4 Stressed

We repeat the same analysis as for stress and also take into account that there are multiple reminders.

my_data.Q.selfc = select_data(c("stressed"))
my_stages = 5:8
my_data.Q.selfc %>% plot_stages(my_stages=my_stages)

Figure 3.48: Embarrassed ratings over time. Lines are mean responses.

the_data = make_the_data(my_data.Q.selfc, phase = "reminder")[, outcome := "selfc"]
attr(the_data,"N")

Table 3.24: Number of particicipants per condition by condition and stage

Drug	state.condition	pre.reminder	post.reminder
oxycodone	control	62	62
oxycodone	stress	63	63
placebo	control	63	63
placebo	stress	63	63

the_data %>% plot_pre_post("reminder") + scale_col_stress

Figure 3.49: Change in ‘stressed’ ratings from pre to post state induction. Each line represents an individual. The stress induction effect is reflected in the slope of the lines.

We can further summarize this information by showing the distribution of the slopes in the previous figure:

plot_stresseff_raw(the_data)

Figure 3.50: Distribution on change in embarrassed ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

plot_stresseff_raw2(the_data)

Figure 3.51: Mean and CIs of change in selfconcious ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Here we run the Bayesian ordinal logistic regression:

fn = "brmsP1/Reminder/stressed_acat.Rdata"
stress_formula = response.c ~ induced.state + Sex + induced.state:Sex + session + (1|participant) + (1|stage)
the_data = the_data %>% 
  coarsen(var = "response", bins = 20, range = c(0,100))
if (file.exists(fn)) {
  load(fn)
} else {
  fit = fit_ordinal_model(the_data, stress_formula, fn, family = acat())
}
#pp_check(fit, ndraws = 2000, type = "bars")

We calculate the contrasts manually from posterior predictions:

tmp = plot_prepost_contrast(fit, the_data,"stressed", my_contrast = "stress")
tmp[[1]]

Figure 3.52: Stress effects in women, men, and sex difference (women-men)

tmp$table

Table 3.25: Modeled post minus pre state induction on stressed ratings by condition (control, stress) and the stress effect by sex.

Sex	control	stress	stress_effect
men	1 [0, 1];	2 [0, 4];	1 [0, 3];
women	1 [0, 3];	3 [0, 6];	2 [0, 3];
total	1 [0, 2];	2 [0, 5];	1 [0, 3];
Sex difference	1 [0, 2];	1 [0, 2];	0 [-0, 1];

post = tmp[[2]]

The estimated effect of stress from the state reminder on feeling stressed is NaN (CI=[NA, NA]; BFNA; P(S>C)NA). The effect is 2 (CI=[0, 3]; BF=41.1; P(S>C)=0.98) in females, 1 (CI=[0, 3]; BF=41.1; P(S>C)=0.98) in males and the sex difference is NaN (CI=[NA, NA]; BFNA; P(S>C)NA).

3.8 Drug effects

3.8.1 High

my_data.Q.high = select_data(c("high"))
my_stages = 1:11
my_data.Q.high %>% 
  .[stage %in% my_stages] %>% 
    .[,.(m = mean(response), se = sd(response)/sqrt(length(unique(participant)))), 
      by = .(stage,Sex,Drug)] %>% 
    .[, `:=`(lower = m-2*se, upper = m+2*se)] %>% 
    ggplot(aes(x = stage, y = m, color = Drug, lty = Sex, fill = Drug)) + 
    geom_line() + 
    geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .2, color = NA) + 
    scale_x_continuous(breaks = my_stages, labels = stage.names[my_stages]) +
  scale_col_drug

Figure 3.53: Feeling high over time. Lines are mean responses.

the_data = make_the_data(my_data.Q.high[stage %in% 4:5],phase = "drug")[, outcome := "high"]
attr(the_data,"N")

Table 3.26: Number of particicipants per condition by condition and stage

Drug	state.condition	pre.drug	post.drug
oxycodone	control	62	62
oxycodone	stress	63	63
placebo	control	63	63
placebo	stress	63	63

the_data %>% plot_pre_post("drug")

Figure 3.54: Change in high from drug ratings from pre to post the first drug administration. Each line represents an individual. The stress induction effect is reflected in the slope of the lines.

We can further summarize this information by showing the distribution of the slopes in the previous figure:

plot_drug_effect_raw(the_data)

Figure 3.55: Distribution on change in high from drug ratings from pre to post state induction, stratified by drug condition (panels) and sex (colors).

Here we run the Bayesian ordinal logistic regression:

fn = "brmsP1/high_acat.Rdata"
the_data = the_data %>%
  .[, Drug := factor(Drug, levels = c("placebo","oxycodone"))] %>% 
  coarsen(var = "response", bins = 20, range = c(0,100))
drug_formula = response.c ~ Drug:stage:Sex + (1 | session) + (1 | participant) 
#drug_formula = response.c ~ Drug:stage:Sex + (1 | session) + (1 + Drug:stage | participant)
if (file.exists(fn)) {
  load(fn)
} else {
  fit = fit_ordinal_model(the_data, drug_formula, fn, family = acat())
}
#pp_check(fit, ndraws = 2000, type = "bars")

Click to show sample diagnostics

Table 3.27: Summary of model paramters

parameter	mean	sd	q5	q95	rhat	ess_bulk	ess_tail
Threshold[1]	1.8	0.73	0.6	3.0	1.006	921	1711
Threshold[2]	0.2	0.76	-1.0	1.4	1.006	1017	1870
Threshold[3]	1.0	0.81	-0.3	2.3	1.005	1108	1948
Threshold[4]	0.5	0.88	-0.9	2.0	1.003	1255	1892
Threshold[5]	-0.2	0.90	-1.6	1.3	1.003	1313	2200
Threshold[6]	0.3	0.88	-1.2	1.7	1.004	1259	2240
Threshold[7]	1.4	0.96	-0.2	3.0	1.004	1516	2189
Threshold[8]	-0.1	1.03	-1.8	1.5	1.002	1578	2244
Threshold[9]	-0.4	0.88	-1.9	1.0	1.003	1319	2340
Threshold[10]	1.3	0.87	-0.1	2.7	1.004	1291	2150
Threshold[11]	-0.2	0.86	-1.6	1.2	1.005	1246	2295
Threshold[12]	1.0	0.84	-0.3	2.4	1.004	1276	2146
Threshold[13]	0.5	0.84	-0.9	1.9	1.004	1195	2080
Threshold[14]	1.2	0.86	-0.2	2.6	1.004	1316	2253
Threshold[15]	0.7	0.87	-0.8	2.1	1.004	1248	2289
Threshold[16]	3.1	1.16	1.3	5.0	1.002	1848	2491
Threshold[17]	-1.4	1.15	-3.3	0.5	1.002	1903	2300
Threshold[18]	1.6	0.93	0.1	3.2	1.002	1526	2282
Threshold[19]	-0.7	0.91	-2.3	0.7	1.005	1351	2256
b_Drugplacebo:stagepre.drug:Sexman	-0.2	0.71	-1.3	1.0	1.006	872	1579
b_Drugoxycodone:stagepre.drug:Sexman	-0.2	0.71	-1.4	0.9	1.006	885	1673
b_Drugplacebo:stagepost.drug:Sexman	-0.2	0.71	-1.3	1.0	1.006	881	1465
b_Drugoxycodone:stagepost.drug:Sexman	0.7	0.71	-0.4	1.9	1.006	900	1637
b_Drugplacebo:stagepre.drug:Sexwoman	-0.4	0.71	-1.6	0.7	1.007	892	1558
b_Drugoxycodone:stagepre.drug:Sexwoman	-0.4	0.72	-1.5	0.8	1.007	895	1377
b_Drugplacebo:stagepost.drug:Sexwoman	-0.2	0.71	-1.3	1.0	1.007	880	1594
b_Drugoxycodone:stagepost.drug:Sexwoman	0.8	0.71	-0.4	1.9	1.007	905	1439
sd_participant__Intercept	0.3	0.05	0.2	0.4	1.000	1591	2265
sd_session__Intercept	0.2	0.22	0.1	0.6	1.001	860	809
model: response.c ~ Drug:stage:Sex + (1 | session) + (1 | participant) b_ = fix effect parameters sd_ standard deviation of hierarchical parameters ess = effective sample size

We calculate the contrasts manually from posterior predictions:

tmp = plot_prepost_contrast(fit, the_data, "high", my_contrast = "drug")
tmp[[1]]

Figure 3.56: Stress effects in women, men, and sex difference (women-men)

tmp$table

Table 3.28: Modeled post minus pre drug administration on high ratings by condition (placebo, oxycodone) and the drug effect by sex.

Sex	placebo	oxycodone	drug_effect
men	1 [-3, 4];	49 [43, 55];	49 [42, 55];
women	2 [0, 5];	58 [53, 63];	56 [50, 62];
total	2 [-1, 4];	54 [50, 58];	52 [48, 57];
Sex difference	2 [-3, 6];	9 [1, 17];	7 [-1, 16];

post = tmp[[2]]

The estimated effect of drug on ratings of being high is NaN (CI=[NA, NA]; BFNA; P(D>P)NA). The effect is 56 (CI=[50, 62]; BF>10^5; P(D>P)>0.999) in females, 49 (CI=[42, 55]; BF>10^5; P(D>P)>0.999) in males and the sex difference is NaN (CI=[NA, NA]; BFNA; P(D>P)NA).

3.9 Drug side effects

3.9.1 Side effect: dizzy

Raw data for ‘dizzy’

Estimated effects of Drug for ‘dizzy’

(#tab:side_effects)Modeled post minus pre drug administration on dizzy ratings by condition (placebo, oxycodone) and the drug effect by sex.

Sex	placebo	oxycodone	drug_effect
men	2 [-1, 5];	34 [28, 40];	32 [26, 39];
women	-1 [-6, 4];	59 [53, 64];	60 [52, 67];
total	0 [-2, 3];	46 [42, 50];	46 [41, 51];
Sex difference	-3 [-9, 3];	24 [16, 32];	27 [18, 37];

3.9.2 Side effect: dry_mouth

Raw data for ‘dry_mouth’

Estimated effects of Drug for ‘dry_mouth’

(#tab:side_effects)Modeled post minus pre drug administration on dry_mouth ratings by condition (placebo, oxycodone) and the drug effect by sex.

Sex	placebo	oxycodone	drug_effect
men	-6 [-12, -1];	-5 [-10, -0];	1 [-6, 8];
women	-2 [-7, 4];	-1 [-6, 5];	1 [-7, 9];
total	-4 [-8, -0];	-3 [-7, 1];	1 [-4, 6];
Sex difference	5 [-3, 12];	4 [-3, 12];	-0 [-11, 10];

3.9.3 Side effect: blunted

Raw data for ‘blunted’

Estimated effects of Drug for ‘blunted’

(#tab:side_effects)Modeled post minus pre drug administration on blunted ratings by condition (placebo, oxycodone) and the drug effect by sex.

Sex	placebo	oxycodone	drug_effect
men	-1 [-7, 4];	34 [27, 40];	36 [27, 44];
women	5 [-0, 9];	49 [43, 55];	45 [37, 52];
total	2 [-2, 5];	42 [37, 46];	40 [34, 45];
Sex difference	6 [-1, 13];	15 [6, 24];	9 [-2, 20];

3.9.4 Side effect: nauseous

Raw data for ‘nauseous’

Estimated effects of Drug for ‘nauseous’

(#tab:side_effects)Modeled post minus pre drug administration on nauseous ratings by condition (placebo, oxycodone) and the drug effect by sex.

Sex	placebo	oxycodone	drug_effect
men	-2 [-5, 1];	4 [1, 8];	6 [2, 11];
women	-2 [-5, 0];	7 [3, 11];	9 [4, 14];
total	-2 [-4, -0];	5 [3, 8];	7 [4, 11];
Sex difference	-0 [-4, 4];	2 [-3, 8];	3 [-4, 9];

3.9.5 Side effect: warm_face

Raw data for ‘warm_face’

Estimated effects of Drug for ‘warm_face’

(#tab:side_effects)Modeled post minus pre drug administration on warm_face ratings by condition (placebo, oxycodone) and the drug effect by sex.

Sex	placebo	oxycodone	drug_effect
men	-6 [-11, -2];	13 [6, 20];	19 [11, 27];
women	-11 [-16, -5];	19 [12, 26];	30 [21, 38];
total	-8 [-12, -5];	16 [11, 21];	24 [18, 30];
Sex difference	-5 [-12, 2];	6 [-4, 16];	11 [-2, 22];

3.9.6 Side effect: not_myself

Raw data for ‘not_myself’

Estimated effects of Drug for ‘not_myself’

(#tab:side_effects)Modeled post minus pre drug administration on not_myself ratings by condition (placebo, oxycodone) and the drug effect by sex.

Sex	placebo	oxycodone	drug_effect
men	-2 [-7, 3];	20 [13, 27];	22 [14, 31];
women	-5 [-11, 1];	34 [27, 41];	40 [30, 49];
total	-4 [-8, 0];	27 [22, 32];	31 [25, 37];
Sex difference	-3 [-11, 5];	14 [4, 24];	17 [4, 30];

3.10 Participants will use vigorous button presses

We display the number of button presses participants made, stratified by drug condition and target dose. Participants with a target dose larger than 50% had to press buttons to achieve a high dose and with a target dose lower than 50% had to press buttons to achieve a low dose.

dt.presses = 
  my_data %>% 
  .[cheater == FALSE] %>% 
  setkeyv(c("participant","session")) %>% 
  .[, last.trial := ifelse(time == max(time),TRUE,FALSE), by = .(session, participant)] %>% 
  .[, n.presses := p_down + p_up] %>% 
  .[, n.presses := cumsum(n.presses), by = .(session, participant)] %>% 
  .[, presses.per.s := n.presses/time] %>% 
  .[last.trial == TRUE] %>% 
  .[!is.na(satisfact)] %>% 
  .[, stage := 0] %>% 
  .[, N.oxy := sum(Drug == "oxycodone"), by = .(participant)] %>% 
  .[N.oxy > 1] %>% 
  .[, .(Sex, participant, session, Drug, 
        state.condition, target, n.presses, presses.per.s)] 


dt.presses[, `target effect` := ifelse(target > 50, "target > 50%","target < 50%")]

hist.presses = 
  dt.presses %>% 
  ggplot(aes(x = presses.per.s, fill = Drug)) + 
  geom_histogram(position = "identity", alpha = .5, bins = 15) + 
  facet_grid(Sex~`target effect`) + 
  xlab("Number of presses second") + 
  ylab("Number of participants") +
  scale_col_drug

hist.presses

dt.presses %>% 
  .[, .(mean = round(mean(presses.per.s),1), sd = round(sd(presses.per.s),1), N = length(unique(participant))),
   by = .(Drug,state.condition)] %>% 
  kable(caption = "Number of button presses per second by drug and stress condition and sex") %>% 
  kable_styling(full_width = FALSE)

Table 3.29: Number of button presses per second by drug and stress condition and sex

Drug	state.condition	mean	sd	N
placebo	stress	5.6	0.5	52
placebo	control	5.5	0.8	50
oxycodone	stress	5.6	0.6	56
oxycodone	control	5.5	0.6	56

dt.presses %>% 
  .[, .(mean = round(mean(n.presses)), sd = round(sd(n.presses))),
   by = .(Drug,state.condition,Sex)] %>% 
  kable(caption = "Number of button presses by drug and stress condition and sex") %>% 
  kable_styling(full_width = FALSE)

Table 3.30: Number of button presses by drug and stress condition and sex

Drug	state.condition	Sex	mean	sd
placebo	stress	women	662	64
placebo	control	women	651	71
oxycodone	stress	women	636	50
oxycodone	control	women	630	84
placebo	control	men	672	123
placebo	stress	men	687	69
oxycodone	control	men	689	51
oxycodone	stress	men	705	75

4 Drug wanting and liking after drug administration

The data for the second drug administration includes only participants who wanted a second dosage larger than zero.

my_data[is.na(satisfact), last.trial := ifelse(time == max(time),TRUE,FALSE), by = .(session, participant)]

my_data_beh = 
  my_data[last.trial == TRUE & is.na(satisfact)] %>% 
  .[, stage := 0] %>% 
  .[, N.oxy := sum(Drug == "oxycodone"), by = .(participant)] %>% 
  .[N.oxy > 1] %>% 
  .[, s.effect_obtained := effect_obtained-mean(effect_obtained), by = .(participant)] %>% 
  .[, s.effect_obtained := s.effect_obtained + mean(effect_obtained)] 
pdata = 
  my_data.Q %>% 
  setkeyv(key(my_data_beh)) %>% 
  .[item_name %in% c("drug_like","drug_again","drug_dislike") & stage %in% c(5,11)] %>% 
  .[my_data, effect_obtained := effect_obtained] %>% 
  .[stage == 5, effect_obtained := 100] %>% 
  .[stage == 11 & effect_obtained == 0, effect_obtained := NA] %>% 
  .[, .(response = mean(response, na.rm = TRUE), sd = sd(response, na.rm = TRUE), N = .N),
    by = .(state.condition,Drug,stage,item_name)] %>% 
  .[, se := sd/sqrt(N)] %>% 
  .[, `:=`(lower=response-2*se,
           upper=response+2*se)] %>% 
  .[, time := stage-4] %>% 
  .[time>1, time := time-4] %>% 
  .[state.condition == "control", time:=time-.05] %>% 
  .[state.condition != "control", time:=time+.05] %>% 
  .[, Stage := factor(stage.names[stage],levels = stage.names)]

pdata %>% 
  ggplot(aes(x = Stage, y = response, ymin=lower, ymax=upper,
             color = Drug, shape = state.condition), group = state.condition) +
  geom_pointrange(size = 1, position = position_dodge(.2)) +
  #geom_vline(xintercept = c(0.75,2.75), lty = 2, alpha = .5) +
  #scale_x_continuous(breaks = c(0.75,2.75), labels = c("dose 1","dose 2")) +  
  ylab("Response") + 
  guides(color=guide_legend(nrow=2,byrow=TRUE),
         shape=guide_legend(nrow=2,byrow=TRUE)) +
  scale_col_drug + 
  xlab("Time in the experiment") +
  facet_wrap(~item_name) +
  theme(legend.position = "bottom") + 
  theme(axis.line = element_line())

Figure 4.1: Drug liking, disliking, and wanting after the administration of the first and second dose.

5 Main Results

5.1 Effect obtained

Prepare data

my_data[is.na(satisfact), last.trial := ifelse(time == max(time),TRUE,FALSE), by = .(session, participant)]

my_data %>% 
  .[, n.presses := p_down + p_up] %>% 
  .[, n.presses := cumsum(n.presses), by = .(session, participant)]

my_data_beh = 
  my_data[last.trial == TRUE & is.na(satisfact)] %>% 
  .[, stage := 0] %>% 
  .[, N.oxy := sum(Drug == "oxycodone"), by = .(participant)] %>% 
  .[N.oxy > 1] %>% 
  .[, s.effect_obtained := effect_obtained-mean(effect_obtained), by = .(participant)] %>% 
  .[, s.effect_obtained := s.effect_obtained + mean(effect_obtained)] 

load("data/dosage_data.Rdata")
feel_tresh = dosage_data[Drug == "placebo", quantile(drug_feel, .95)]
dosage_data = 
  dosage_data %>% 
  .[, session := as.integer(session)] %>% 
  setkeyv(.,key(my_data_beh)) %>% 
  .[my_data_beh, target := target] %>% 
  .[, drug_2_refused := 
      ifelse(is.na(drug_2_admin) & 
               drug_2_dose == 0 & 
               target > 0 & 
               drug_feel < feel_tresh, TRUE, FALSE)]

dosage.factor = 
  median(dosage_data$weight_kg/dosage_data$drug_1_dose,
       na.rm = TRUE)

mg_per_kg = round(mean(dosage_data$drug_1_dose/dosage_data$weight_kg),3)

m.format = function(x) round(mean(x,na.rm = TRUE),2)

5.1.1 Drug dosages

The first oxycodone dose was 0.043 mg per kg body weight, and an average dose of 3.12 mg. The dose for the second drug administration was determined by the result of the self administration task, whereby participant who had worked to obtain more than 100% of the effect after the first dose still received 100%. In the oxycodone condition after stress induction participants received on average 2.34 mg at the second administration (working for on average 85.02%), whereas they received after control induction on average 2.16 mg (working for on average 77.68%). Considering only those who wanted a non-zero dose, participants received on average 2.88 mg n the oxycodone condition after stress induction at the second administration (working for on average 97.4%), whereas they received after control induction on average 2.62 mg (working for on average 92.37%).

3 participants, who had worked to obtain on average 3.75% of the first effect, decided that they did not want the second drug administration. 13 participants worked to not obtain any drug at the second administration.

Drug effects split by Sex:

my_data_beh.stat = 
  my_data_beh %>% 
  .[, .(effect_obtained = mean(effect_obtained),
        se = sd(effect_obtained)/sqrt(.N),
        N = .N),
    by = .(state.condition,Drug,stage,Sex)] %>% 
  .[, `:=`(lower = effect_obtained-2*se,upper = effect_obtained+2*se)]

my_data_beh.stat %>% 
  ggplot(aes(x = state.condition, y = effect_obtained, color = Drug, group = Drug)) + 
  geom_quasirandom(data = my_data_beh, alpha = .25) + 
  geom_line(position = position_dodge(.5)) +
  geom_point(position = position_dodge(.5), size = 4) + 
  geom_linerange(aes(ymin = lower, ymax = upper), position = position_dodge(.50)) + 
  facet_grid(Drug~Sex) +
  scale_col_drug + 
  ylim(0,126) +
  theme(panel.spacing = unit(2, "lines")) + 
  ylab("Dose obtained for the 2nd drug administration")

Figure 5.1: Obtained dosage by drug, state induction and sex.

my_data_beh %>% 
  .[, .(participant,Sex, state.condition, effect_obtained, Drug)] %>% 
  dcast(participant + Drug + Sex ~state.condition, value.var = "effect_obtained") %>% 
  .[, stress_effect := stress-control] %>% 
  ggplot(aes(x = Sex, y = stress_effect, color = Sex, fill = Sex)) + 
  scale_col_sex + 
  stat_summary(fun.y = mean,alpha = .75,
               geom = "bar", color = NA) + 
  stat_summary(fun.data = mean_cl_normal, color = "black",
               geom = "errorbar", width = .25) + 
  facet_wrap(~Drug) + 
  geom_quasirandom(shape = 1, color = "black") + 
  theme(legend.position = "right")

Figure 5.2: Effect of stress on obtained effect by sex and drug in raw data.

my_data_beh %>% 
  .[, .(participant,Sex, state.condition, effect_obtained, Drug)] %>% 
  dcast(participant + Drug + Sex ~state.condition, value.var = "effect_obtained") %>% 
  .[, stress_effect := stress-control] %>% 
  dcast(participant + Sex ~ Drug, value.var = "stress_effect") %>% 
  .[, stress_effect_x_drug := oxycodone - placebo] %>% 
  ggplot(aes(x = Sex, y = stress_effect_x_drug, color = Sex, fill = Sex)) + 
  scale_col_sex + 
  stat_summary(fun.y = mean,alpha = .75,
               geom = "bar", color = NA) + 
  stat_summary(fun.data = mean_cl_normal, color = "black",
               geom = "errorbar", width = .25) + 
  geom_quasirandom(shape = 1, color = "black") + 
  theme(legend.position = "right")

Figure 5.3: Interaction of the effct of stress on obtained effect with drug by sex. Larger values indicate greater drug-induced wanting increase in the ocyxodone condition.

Raw data: Select only the oxycodone condition (According to pre-registered plan):

my_data_beh.all = 
  my_data_beh[, .(Sex, participant, session,
                  Drug, state.condition, effect_obtained, target,
                  mean.delta_effect_obtained,
                  first.time.target,n.presses)] 
my_data_beh = 
  copy(my_data_beh.all) %>% 
  .[Drug == "oxycodone"]

To deal with the multi-modal nature of the data, we use an ordered logistic regression, for which we coarsen the data in 25 equally sized ordered categories.

n_cat = 25
my_data_beh = 
  coarsen(my_data_beh,
          var = "effect_obtained",
          bins = n_cat,
          range = c(0,125))

5.1.2 Linear regression without sex effects

We implement a hierarchical regression model where random effects account for repeated measurements. Even though the raw data show that the effect of stress is stronger in men, we first analyse the data without modelling sex effects. The regression model uses fixed effects parameters to model the effects of stress manipulation and random effects to model varying intercepts and stress effects across participants and effects of session number.

my_data_beh[, session := factor(session)]

We estimate a standard hierarchical linear regression with Gaussian errors.

ffit = lme4::lmer(
  effect_obtained ~ state.condition + (1 | participant) + (1|session),
  data = my_data_beh)

To check the model, we first investigate model predictions. The following figure shows that the distribution of model predictions deviates strongly from the observed values and also include impossible negative values.

error_sd = 
  VarCorr(ffit) %>% attr("sc")
p1 = 
  data.table(yhat = predict(ffit)) %>% 
  .[, lower := yhat - error_sd] %>% 
  .[, upper := yhat + error_sd] %>% 
  .[, state.condition := my_data_beh$state.condition] %>% 
  .[order(yhat)] %>% 
  .[, j := 1:.N] %>% 
  ggplot(aes(x = j, y = yhat, ymin = lower, ymax = upper,
             color = state.condition, fill = state.condition)) + 
   geom_hline(yintercept = 0, color = "red") +
  geom_linerange(alpha = .5) +
  geom_point(size = .4) +
  xlab("") +
  ylab("fitted + error") + 
  scale_col_stress
p2 = 
  rbind(my_data_beh[, .(effect_obtained)][, data := "observed"],
           data.table(effect_obtained = predict(ffit), data = "predicted")) %>% 
  ggplot(aes(x = effect_obtained, fill = data)) + 
  geom_histogram(position = "identity", alpha = .5)
  
(p1 | p2) + plot_layout(widths = c(2,1))

Figure 5.4: Predictions of the linear model +/- 1sd of the normal error distribution.

Consistently, the following figure shows that the residuals are correlated with the outcome variable and not normally distributed.

res = data.table(residuals = residuals(ffit), fitted = predict(ffit))
p1 = res %>% 
  ggplot(aes(x = fitted, y = residuals)) + 
  geom_point() + 
  geom_smooth()
p2 = 
  res %>% 
  ggplot(aes(x = residuals)) + 
  geom_histogram(aes(y =..density..),
                 breaks = seq(-60, 60, by = 5), 
                 colour = "black", 
                 fill = "white") +
  stat_function(fun = dnorm,
                args = list(mean = mean(res$residuals),
                            sd = sd(res$residuals)))

p1 | p2

Figure 5.5: Residuals from a linear regression model. Left: Predicted values and residuals are not independent. Right: Residuals are not normally distributed.

rm(error_sd,ffit,p1,p2,res)

5.1.3 Ordinal logistic regression (OLR) without sex effects

As for the questionnaire data a linear regression model with Gaussian errors is not a good model for the obtained dosage data. We therefore use ordinal regression models for coarsened data. Specifically, we coarsen the obtained dosages data into 20 equally sized bins and perform Bayesian ordinal logistic regressions on the coarsened data. To report results, we back-calculate effects (including credible intervals) of independent variables on the original 0-125% scale.

We implement a hierarchical ordinal regression model where the cumulative logistic models accounts for the multi-modal distribution of the obtained drug effects and random effects account for repeated measurements. Even though the raw data show that the effect of stress is stronger in men, we first analyse the data without modelling sex effects. The regression model uses fixed effects parameters to model the effects of stress manipulation and random effects to model varying intercepts and stress effects across participants and effects of session number.

analysis_formula = 
  formula(effect_obtained.c ~ state.condition + (1 + state.condition | participant) + (1|session))

We use brms’ default priors for standard deviations of random effects (student-t with 3 df and a sd of 2.5) and a weakly informative normal prior for main effects (normal distribution with mean 0 and standard deviation 2). Priors for the thresholds are specified as normal distributions with the logistic quantile function at cumulative probabilities of response categories as means and standard deviations of 1.

fn = "brmsP1/obtained/nosex.Rdata"
if (file.exists(fn)) {
  load(fn)
} else {
  my_prior = add_cumthresh_prior(
    prior(normal(0,2), class = b),
    my_data_beh$effect_obtained.c)
  fit = brm(
    analysis_formula,
    family = cumulative(),
    data = my_data_beh,
    prior = my_prior,
    backend = "cmdstanr",
    iter = 2000)
  save(fit,file = fn)
}

Click to show sample diagnostics

Table 5.1: Summary of model paramters

parameter	mean	sd	q5	q95	rhat	ess_bulk	ess_tail
b_Intercept[1]	-3.0	0.38	-3.6	-2.4	1.002	4455	2550
b_Intercept[2]	-2.6	0.35	-3.2	-2.0	1.000	5015	3363
b_Intercept[3]	-2.3	0.32	-2.9	-1.8	1.000	4951	2721
b_Intercept[4]	-2.2	0.32	-2.7	-1.7	1.000	5225	2861
b_Intercept[5]	-2.1	0.31	-2.6	-1.6	1.000	5232	2949
b_Intercept[6]	-1.9	0.30	-2.4	-1.5	1.001	5444	2772
b_Intercept[7]	-1.8	0.30	-2.3	-1.3	1.001	5484	2800
b_Intercept[8]	-1.5	0.29	-2.0	-1.1	1.000	5556	3330
b_Intercept[9]	-1.4	0.29	-1.9	-1.0	1.000	5490	3226
b_Intercept[10]	-1.3	0.29	-1.8	-0.9	1.000	5770	3017
b_Intercept[11]	-1.2	0.29	-1.6	-0.7	1.000	5282	2919
b_Intercept[12]	-1.0	0.29	-1.5	-0.5	1.000	5522	3078
b_Intercept[13]	-0.8	0.29	-1.3	-0.3	1.001	5332	3133
b_Intercept[14]	-0.6	0.29	-1.1	-0.2	1.000	5555	3412
b_Intercept[15]	-0.4	0.30	-0.9	0.1	1.001	5290	3649
b_Intercept[16]	-0.1	0.30	-0.6	0.4	1.001	4803	3383
b_Intercept[17]	0.7	0.33	0.2	1.2	1.000	3862	3159
b_Intercept[18]	1.1	0.35	0.5	1.7	1.000	3468	2951
b_Intercept[19]	1.4	0.36	0.8	2.0	1.000	3172	2536
b_Intercept[20]	1.5	0.38	0.9	2.2	1.000	2955	2102
b_state.conditionstress	0.5	0.36	-0.1	1.1	1.001	4779	2654
sd_participant__Intercept	2.1	0.46	1.4	2.9	1.003	930	1523
sd_participant__state.conditionstress	0.7	0.56	0.1	1.8	1.003	793	1007
sd_session__Intercept	0.8	0.53	0.1	1.7	1.002	1410	1672
model: effect_obtained.c ~ state.condition + (1 + state.condition | participant) + (1 | session) b_ = fix effect parameters sd_ standard deviation of hierarchical parameters ess = effective sample size

The next figure shows that the model describes the observed data well.

pp_check(fit, ndraws = 2000, type = "bars")

Figure 5.6: Posterior predictive check for obtained effects.

5.1.4 OLR without sex effects, adjacent category model with category specific effects

We repeat the same analysis as before, except that we use an adjacent category ordered logistic regression with category specific effects. Differently than the previous analysis, this model does not assume proportional odds. The advantage of this model is that we can fit the data better and have a higher power to detect effects. The disadvantage is that the model estimates for each category transition separate odds ratios. Therefore, it is not meaningful to report one odds ratio as a model parameter.

Here we estimate the model.

analysis_formula = 
  formula(effect_obtained.c ~ cs(state.condition) + (1 + state.condition | participant) + (1|session))
prior_sds = c(.5, 1, 5, 2)
make_prior_b = function(prior_sd) {
  return(
    eval(parse(text = paste("prior(normal(0,",prior_sd,"), class = b)")))
  )
}
for (prior_sd in prior_sds) {
  fn = paste0("brmsP1/obtained/nosex_acat",prior_sd,".Rdata")
  if (file.exists(fn)) {
    load(fn)
  } else {
    my_prior = 
      get_acat_thresh_priors(
        my_data_beh$effect_obtained.c,
        make_prior_b(prior_sd) + 
          prior(normal(0,2), class = "sd"))
    
    fit = brm(
      analysis_formula,
      family = acat(),
      data = my_data_beh,
      prior = my_prior,
      backend = "cmdstanr",
      iter = 2000)
    save(fit,file = fn)
  }
}

Click to show sample diagnostics

Table 5.2: Summary of model paramters

parameter	mean	sd	q5	q95	rhat	ess_bulk	ess_tail
Threshold[1]	0.7	1.03	-1.0	2.4	1.001	1834	2639
Threshold[2]	-2.0	1.11	-3.9	-0.2	1.000	2012	2413
Threshold[3]	0.1	1.20	-1.8	2.2	1.001	3144	2679
Threshold[4]	-1.3	1.34	-3.6	0.9	1.003	3554	2433
Threshold[5]	-0.9	1.33	-3.1	1.3	1.002	2963	2455
Threshold[6]	-2.0	1.25	-4.1	-0.1	1.000	3149	2710
Threshold[7]	-1.8	1.05	-3.6	-0.2	1.002	2764	2499
Threshold[8]	0.5	1.10	-1.2	2.4	1.000	3333	2646
Threshold[9]	-0.7	1.24	-2.8	1.3	1.002	3640	2867
Threshold[10]	-1.7	1.11	-3.6	0.0	1.001	3569	2704
Threshold[11]	-0.2	0.96	-1.8	1.4	1.001	3173	2532
Threshold[12]	-0.7	0.95	-2.2	0.9	1.000	3067	2509
Threshold[13]	-0.5	0.90	-2.0	1.0	1.001	2997	2803
Threshold[14]	-0.5	0.83	-1.9	0.9	1.001	3224	2590
Threshold[15]	0.8	0.82	-0.6	2.1	1.001	2912	2238
Threshold[16]	-1.1	0.79	-2.4	0.2	1.001	2140	2377
Threshold[17]	1.2	0.71	0.1	2.4	1.002	1737	2377
Threshold[18]	1.3	0.85	0.0	2.8	1.001	2039	2048
Threshold[19]	2.1	1.04	0.5	3.9	1.001	2007	2476
Threshold[20]	-1.4	1.02	-3.1	0.3	1.003	1351	2180
bcs_state.conditionstress[1]	2.2	1.22	0.2	4.3	1.002	2657	2932
bcs_state.conditionstress[2]	-0.5	1.33	-2.7	1.7	1.000	3301	2684
bcs_state.conditionstress[3]	1.1	1.45	-1.3	3.4	1.000	3443	2927
bcs_state.conditionstress[4]	-0.3	1.55	-2.9	2.2	1.001	3714	2947
bcs_state.conditionstress[5]	1.4	1.56	-1.1	4.0	1.000	3775	2714
bcs_state.conditionstress[6]	0.1	1.43	-2.3	2.4	1.001	3522	2681
bcs_state.conditionstress[7]	0.2	1.29	-1.9	2.3	1.001	3582	2999
bcs_state.conditionstress[8]	-0.8	1.48	-3.3	1.5	1.001	4293	2917
bcs_state.conditionstress[9]	0.3	1.59	-2.3	2.9	1.000	4269	2628
bcs_state.conditionstress[10]	1.3	1.49	-1.1	3.8	1.001	4254	3015
bcs_state.conditionstress[11]	1.2	1.27	-0.8	3.3	1.002	3915	3206
bcs_state.conditionstress[12]	-0.5	1.26	-2.6	1.6	1.000	3674	3088
bcs_state.conditionstress[13]	-0.7	1.34	-2.9	1.5	1.000	4276	3131
bcs_state.conditionstress[14]	1.3	1.34	-0.9	3.5	1.001	3599	3095
bcs_state.conditionstress[15]	2.4	1.11	0.6	4.3	1.001	3367	2842
bcs_state.conditionstress[16]	-0.2	0.89	-1.7	1.2	1.000	2910	2595
bcs_state.conditionstress[17]	-0.5	0.87	-1.9	1.0	1.001	3826	3131
bcs_state.conditionstress[18]	-0.5	1.12	-2.3	1.3	1.001	3200	2915
bcs_state.conditionstress[19]	0.5	1.29	-1.5	2.7	1.001	3402	2874
bcs_state.conditionstress[20]	0.3	1.19	-1.6	2.3	1.001	3009	2847
sd_participant__Intercept	1.6	0.55	0.8	2.6	1.007	528	763
sd_participant__state.conditionstress	1.8	0.59	0.9	2.9	1.004	678	1423
sd_session__Intercept	0.6	0.47	0.1	1.6	1.014	635	930
model: effect_obtained.c ~ cs(state.condition) + (1 + state.condition | participant) + (1 | session) b_ = fix effect parameters sd_ standard deviation of hierarchical parameters ess = effective sample size

Next, we calculate the effects on the original 0-125 scale from the fitted model. This approach removes the effect of session by calculating model-predicted outcomes for the two stress conditions in under the assumption that each condition could be in session 1, 2, 3, or 4.

my_contrasts = plot_results_origscale2(fit)

my_contrasts$plot

Figure 5.7: Estimated mean drug wanting in the 2 state induction conditions (left) and effect of stress, calculate as drug wanting in the stress condition minus drug wanting in the conrol condition (left)

my_caption = "Model-estimated effect obtained by condition"
my_contrasts$simple.stats %>% 
  kable(caption = my_caption) %>% 
  kable_styling(full_width = FALSE) 

Table 5.3: Model-estimated effect obtained by condition

state.condition	m (CrI)
stress	89 (86, 92)
control	84 (80, 87)

In the stress condition participants worked on average for 5 (CI = [1, 10], P_E>0 = 0.99) percentage points higher oxycodone dosage.

5.1.5 Prior sensitivity analysis

To investigate the effect of the model priors on the results we compare the effects calculated with the posterior distribution of the parameters with effects calculate from the prior distribution of the parameters.

First, we generate samples from the prior distribution only and we calculate the same contrasts that we previously calculated with the posterior distributions, this time using the prior distribution.

psa_contrasts = c()
for (prior_sd in prior_sds) {
  fn = paste0("brmsP1/obtained/nosex_acat",prior_sd,"_psa.Rdata")
  if (file.exists(fn)) {
    load(fn)
  } else {
    my_prior = 
      eval(parse(text = paste("prior(normal(0,",prior_sd,"), class = b)")))+ 
      prior(normal(0,1), class = Intercept)
    prior_samples = brm(
      analysis_formula,
      sample_prior = "only", # only prior samples
      family = acat(),
      data = my_data_beh,
      prior = my_prior,
      backend = "cmdstanr",
      iter = 5000)
    save(prior_samples,file = fn)
  }
  # calculate contrasts
  my_contrasts_prior = plot_results_origscale2(prior_samples)
  load(paste0("brmsP1/obtained/nosex_acat",prior_sd,".Rdata"))
  my_contrasts = plot_results_origscale2(fit)
  
  # collect all samples
  psa_contrasts = rbind(
    psa_contrasts,
    rbind(
      my_contrasts$dt.samples.d[, .(value)][, from := "posterior"][, prior_sd := prior_sd],
      my_contrasts_prior$dt.samples.d[, .(value)][, from := "prior"][, prior_sd := prior_sd]
    )
  )
}

And we compare stress effects from posterior and prior under the chosen prior of \(\textrm{Normal}(0,2)\) for the effect of state condition.

psa_contrasts %>% 
  .[, from := factor(from, levels = c("prior","posterior"))] %>% 
  .[prior_sd == 2] %>% 
  ggplot(aes(x = value, color = from)) + 
  geom_density() + 
  xlab("calculated effect") +
  ylab("") +
  theme(legend.position = "right") +
  gg_no_y_axis

Figure 5.8: Effect distribution under posterior samples and prior prior samples.

The figure above illustrates that the prior distribution for the effect is a shrinkage prior that favors null effects over alternative effects. In comparison, the posterior distribution of the effect is narrower and has a non-zero mean. Generally, the employed shrinkage prior means that our estimate of the stress effect is conservative.

Here are calculated effects given the prior and posterior distribution for alternative standard deviations for the prior on the state-condition effect:

psa_contrasts %>% 
  .[, `prior sd` := factor(prior_sd, levels = c("2","0.5","1","5"))] %>% 
  ggplot(aes(x = value, color = from, lty = `prior sd`)) + 
  geom_density() + 
  ylab("calculated effect") +
  xlab("") +
  facet_wrap(~from, scales = "free", nrow = 2) +
  theme(legend.position = "right") +
  gg_no_y_axis

Figure 5.9: Effect distribution under posterior samples and prior prior samples.

psa_contrasts %>% 
  .[, .(`mean (CI)` = get_mci(value,digits = 1)), by = .(from,prior_sd)] %>% 
  dcast(prior_sd ~ from) %>% 
  kable(caption = "Effects of stress condition caculated from prior and posterior distribution.") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic()

Table 5.4: Effects of stress condition caculated from prior and posterior distribution.

prior_sd	prior	posterior
0.5	-0.3 [-21.3, 19.9]; 0.4849	4.2 [-0.2, 8.7]; 0.97
1.0	-0.3 [-21.8, 20.0]; 0.4896	4.7 [0.1, 9.2]; 0.9777
2.0	0.2 [-23.0, 24.1]; 0.5097	5.1 [0.8, 9.5]; 0.991
5.0	0.8 [-34.5, 35.4]; 0.521	5.7 [1.5, 10.2]; 0.99475

Priors with smaller standard deviations for the effect of stress lead to posterior distributions for the effect of stress with somewhat smaller means and standard deviations, whereas priors with larger standard deviations for the effect of stress would lead to posterior distributions for the effect of stress with somewhat larger mean and standard deviations. In sum, our results show that we have employed shrinkage priors that result in a conservative estimate of the effect of stress.

5.1.6 OLR with sex effects

Because the initial data visualisation also shows sex effects, we estimate sex specific main effects and effects of stress in the next model.

This analysis uses an adjacent category model with category-specific effects, which relaxes the assumption of proportional odds. The basic analysis model is to estimate the the obtained effect dependent on state-condition and Sex, while accounting for repeated measurement within individuals and for session effects with random effects.

In particular we use the following model:

analysis_formula = 
  formula(effect_obtained.c ~ 
            # category specific effects of state condition and sex
            cs(state.condition*Sex) +
            # random effects for participants
            (1 + state.condition | participant) + 
            # random effects for session
            (1|session))

for (prior_sd in prior_sds) {
  fn = paste0("brmsP1/obtained/acat",prior_sd,".Rdata")
  if (file.exists(fn)) {
    load(fn)
  } else {
    my_prior = 
      get_acat_thresh_priors(
        my_data_beh$effect_obtained.c,
        make_prior_b(prior_sd) + 
          prior(normal(0,2), class = "sd"))
    control = list(adapt_delta = .90)
    fit = brm(
      analysis_formula,
      family = acat(),
      data = my_data_beh,
      prior = my_prior,
      backend = "cmdstanr",
      iter = 2000)
    save(fit,file = fn)
    prior_samples = update(fit,sample_prior = "only")
    save(prior_samples,file = gsub(".Rdata","_psa.Rdata",fn))
  }
}
load(paste0("brmsP1/obtained/acat",2,".Rdata"))

Click to show sample diagnostics

Table 5.5: Summary of model paramters

parameter	mean	sd	q5	q95	rhat	ess_bulk	ess_tail
Threshold[1]	0.9	1.17	-1.0	2.9	1.001	2812	1945
Threshold[2]	-2.4	1.26	-4.5	-0.4	1.003	2357	2792
Threshold[3]	0.7	1.35	-1.5	2.9	1.000	2954	2762
Threshold[4]	-1.0	1.45	-3.4	1.3	1.000	2901	3177
Threshold[5]	-0.9	1.49	-3.4	1.6	1.001	3218	2476
Threshold[6]	-2.2	1.40	-4.5	0.1	1.000	3216	2782
Threshold[7]	-2.4	1.28	-4.6	-0.4	1.000	3112	3274
Threshold[8]	0.9	1.29	-1.2	3.0	1.001	3164	2542
Threshold[9]	-0.6	1.36	-2.8	1.7	1.001	2556	2173
Threshold[10]	-2.0	1.29	-4.1	0.1	1.002	2900	2677
Threshold[11]	-0.2	1.19	-2.1	1.8	1.001	2680	2677
Threshold[12]	-1.2	1.19	-3.2	0.7	1.001	2702	2501
Threshold[13]	0.1	1.11	-1.7	1.9	1.002	2524	2963
Threshold[14]	-1.0	1.08	-2.8	0.8	1.003	2435	2681
Threshold[15]	0.8	1.04	-0.9	2.5	1.003	2244	2611
Threshold[16]	-1.4	0.97	-3.1	0.1	1.001	2010	2748
Threshold[17]	0.9	0.85	-0.5	2.3	1.004	1820	2086
Threshold[18]	1.2	1.01	-0.4	2.9	1.001	1875	2556
Threshold[19]	2.2	1.17	0.4	4.2	1.001	2026	2038
Threshold[20]	-2.2	1.15	-4.1	-0.3	1.003	1481	2029
bcs_state.conditionstress[1]	1.3	1.38	-1.0	3.6	1.001	2661	2492
bcs_state.conditionstress[2]	-0.1	1.47	-2.5	2.4	1.000	3571	2944
bcs_state.conditionstress[3]	0.2	1.59	-2.5	2.8	1.000	3419	2726
bcs_state.conditionstress[4]	-0.8	1.67	-3.6	1.9	1.001	3613	2578
bcs_state.conditionstress[5]	0.8	1.60	-1.9	3.4	1.001	3534	2640
bcs_state.conditionstress[6]	-0.2	1.63	-3.0	2.5	1.001	3833	2801
bcs_state.conditionstress[7]	0.2	1.56	-2.4	2.7	1.002	3205	2992
bcs_state.conditionstress[8]	-0.7	1.65	-3.5	2.0	1.000	4074	3000
bcs_state.conditionstress[9]	0.2	1.67	-2.5	3.0	1.004	4024	2986
bcs_state.conditionstress[10]	1.1	1.69	-1.6	3.9	1.000	4231	2695
bcs_state.conditionstress[11]	1.1	1.57	-1.4	3.7	1.000	3540	2658
bcs_state.conditionstress[12]	0.3	1.51	-2.1	2.8	1.000	3193	2998
bcs_state.conditionstress[13]	0.0	1.55	-2.5	2.5	1.001	3452	2586
bcs_state.conditionstress[14]	1.8	1.54	-0.8	4.3	1.002	3642	3019
bcs_state.conditionstress[15]	2.9	1.34	0.7	5.1	1.001	3345	2610
bcs_state.conditionstress[16]	0.9	1.10	-0.9	2.7	1.002	2714	2697
bcs_state.conditionstress[17]	-0.2	1.01	-1.9	1.4	1.001	2480	2407
bcs_state.conditionstress[18]	-0.8	1.24	-2.8	1.3	1.001	2381	2393
bcs_state.conditionstress[19]	0.5	1.36	-1.7	2.7	1.000	3201	2689
bcs_state.conditionstress[20]	0.8	1.31	-1.4	2.9	1.002	2719	2803
bcs_Sexwoman[1]	0.5	1.28	-1.7	2.5	1.001	2753	2532
bcs_Sexwoman[2]	-0.2	1.39	-2.5	2.1	1.000	2940	3041
bcs_Sexwoman[3]	1.2	1.51	-1.3	3.7	1.002	3030	2175
bcs_Sexwoman[4]	0.9	1.53	-1.6	3.4	1.002	2870	2645
bcs_Sexwoman[5]	0.0	1.58	-2.6	2.5	1.002	3070	2666
bcs_Sexwoman[6]	0.1	1.49	-2.4	2.5	1.000	3428	3054
bcs_Sexwoman[7]	-1.0	1.41	-3.3	1.3	1.001	2910	2964
bcs_Sexwoman[8]	1.1	1.51	-1.4	3.6	1.001	3183	2893
bcs_Sexwoman[9]	0.3	1.51	-2.2	2.8	1.001	2925	2835
bcs_Sexwoman[10]	-0.5	1.46	-2.9	1.9	1.002	2482	2455
bcs_Sexwoman[11]	0.0	1.34	-2.2	2.2	1.002	2953	2868
bcs_Sexwoman[12]	-1.2	1.36	-3.4	1.0	1.001	2923	2849
bcs_Sexwoman[13]	1.2	1.29	-0.8	3.4	1.001	2524	2911
bcs_Sexwoman[14]	-1.1	1.22	-3.1	0.9	1.003	2678	2612
bcs_Sexwoman[15]	-0.3	1.17	-2.1	1.7	1.001	2503	3028
bcs_Sexwoman[16]	-0.8	1.11	-2.6	1.0	1.002	2425	2668
bcs_Sexwoman[17]	-1.1	1.03	-2.9	0.6	1.002	2200	2769
bcs_Sexwoman[18]	-1.0	1.18	-3.0	0.9	1.001	2638	2895
bcs_Sexwoman[19]	-0.7	1.29	-2.8	1.4	1.001	2069	2778
bcs_Sexwoman[20]	-2.7	1.25	-4.8	-0.7	1.001	2394	2678
bcs_state.conditionstress:Sexwoman[1]	1.2	1.44	-1.1	3.6	1.001	2818	2857
bcs_state.conditionstress:Sexwoman[2]	-1.5	1.61	-4.2	1.0	1.000	3049	2846
bcs_state.conditionstress:Sexwoman[3]	1.2	1.66	-1.6	3.9	1.000	4037	3012
bcs_state.conditionstress:Sexwoman[4]	-0.2	1.66	-3.0	2.4	1.000	3621	2821
bcs_state.conditionstress:Sexwoman[5]	1.1	1.71	-1.7	3.9	1.001	3247	2751
bcs_state.conditionstress:Sexwoman[6]	-0.1	1.64	-2.7	2.6	1.002	3353	2489
bcs_state.conditionstress:Sexwoman[7]	0.2	1.60	-2.4	2.8	1.000	3382	2712
bcs_state.conditionstress:Sexwoman[8]	-1.3	1.69	-4.1	1.4	1.001	3683	2905
bcs_state.conditionstress:Sexwoman[9]	-0.4	1.72	-3.3	2.4	1.002	3972	2711
bcs_state.conditionstress:Sexwoman[10]	0.5	1.72	-2.4	3.3	1.001	4229	3165
bcs_state.conditionstress:Sexwoman[11]	0.2	1.60	-2.5	2.8	1.003	3218	2645
bcs_state.conditionstress:Sexwoman[12]	-1.0	1.53	-3.5	1.5	1.000	3321	2972
bcs_state.conditionstress:Sexwoman[13]	-1.9	1.60	-4.5	0.7	0.999	3896	2985
bcs_state.conditionstress:Sexwoman[14]	-0.3	1.62	-3.0	2.4	1.000	4040	3062
bcs_state.conditionstress:Sexwoman[15]	-0.4	1.48	-2.9	2.0	1.002	3331	3084
bcs_state.conditionstress:Sexwoman[16]	-2.1	1.28	-4.2	0.0	1.000	2797	2890
bcs_state.conditionstress:Sexwoman[17]	-1.3	1.34	-3.5	0.8	1.002	3194	2663
bcs_state.conditionstress:Sexwoman[18]	1.3	1.56	-1.3	3.9	1.000	3340	2396
bcs_state.conditionstress:Sexwoman[19]	0.0	1.58	-2.5	2.7	1.001	3533	2876
bcs_state.conditionstress:Sexwoman[20]	-1.0	1.46	-3.4	1.4	1.001	3604	3202
sd_participant__Intercept	1.9	0.59	1.1	3.0	1.005	406	1003
sd_participant__state.conditionstress	2.0	0.67	1.1	3.2	1.009	413	719
sd_session__Intercept	0.9	0.60	0.2	2.0	1.003	486	725
model: effect_obtained.c ~ cs(state.condition * Sex) + (1 + state.condition | participant) + (1 | session) b_ = fix effect parameters sd_ standard deviation of hierarchical parameters ess = effective sample size

To communicate the results we calculate conditional effects, which show by how many percentage points the obtained effects increased on average in the stress condition. See the function make_plot_samples in the utils.R file for details.

plot.stats_table_obtained = 
  plot_results_origscale(
    fit,
    my_data_beh,
    outcome.var = "effect_obtained.c")

p1 = plot.stats_table_obtained$p1
yrange = layer_scales(p1)$y$range$range
p1 = p1 + scale_y_continuous(limits = yrange) +
  guides(y = guide_axis_truncated(
    trunc_lower = yrange[1]+diff(yrange)/500)) +
  coord_cartesian(ylim = c(67,108.5)) +
  xlab("") +
  theme(plot.title = element_text(hjust = -.05))

p2 = plot.stats_table_obtained$p2

px = 
  (p1 | p2) + 
  plot_layout(widths = c(4,1),) +
  plot_annotation(tag_levels = 'A')
ggsave(px,
        filename = "figures/Effect_Obtained.png", 
        width = 20,height = 15, units = "cm")
ggsave(px,
        filename = "figures/Effect_Obtained.pdf", 
        width = 20,height = 15, units = "cm", device = cairo_pdf)

“Effect of stress on effect obtained in the self administration task”

my_caption = "Model-estimated effect obtained by sex and state condition."
my_footnote = "Numbers are means and CIs: Upper and lower bound of 95% credible intervals"
plot.stats_table_obtained$simple.stats %>% 
  kable(digits = 2, caption = my_caption) %>% 
  kable_styling(full_width = FALSE) %>%
  kable_classic() %>% 
  add_footnote(my_footnote, notation="none")

Table 5.6: Model-estimated effect obtained by sex and state condition.

Sex	control	stress
men	91 (86, 96)	104 (99, 108)
total	84 (81, 87)	89 (86, 92)
women	77 (72, 81)	74 (69, 78)
Numbers are means and CIs: Upper and lower bound of 95% credible intervals

plot.stats_table_obtained$tbl_me_sex

Table 5.7: Effect of sex, calculated as effect obtained in men - women.

state condition	effect_woman
control	14 [8, 20]; >.999
stress	30 [24, 37]; >.999
overall	22 [18, 27]; >.999
Numbers are means and CIs: Upper and lower bound of 95% credible intervals

my_caption = "Effect of stress on effect obtained by sex, over all, and sex difference"
my_footnote = "CI: Upper and lower bound of 95% credible intervals. P > 0: Posterior probability that the parameter or contrast is larger than 0."
plot.stats_table_obtained$stats[,.(Sex, mean, CI, `P > 0`)] %>% 
  .[, Sex := factor(Sex, levels = c("total","women","men","sex diff: m-f"))] %>% 
  setkeyv("Sex") %>% 
  kable(digits = 2, caption = my_caption) %>% 
  kable_styling(full_width = FALSE) %>%
  kable_classic() %>% 
  add_footnote(my_footnote, notation="none")

Table 5.8: Effect of stress on effect obtained by sex, over all, and sex difference

Sex	mean	CI	P > 0
total	4.91	0.56, 9.22	0.99
women	-3.07	-9.1, 3.33	0.16
men	12.89	6.72, 18.66	1.00
sex diff: m-f	15.96	7.02, 24.49	1.00
CI: Upper and lower bound of 95% credible intervals. P > 0: Posterior probability that the parameter or contrast is larger than 0.

PS = plot.stats_table_obtained$stats

We used hierarchical ordinal logistic regression models to estimate the effect of stress on the drug wanting, which was measured as the percent of the drug effect from the first drug administration participants worked for in the self administration task. To better understand the magnitude of these effects, we calculated marginal effects of stress on drug wanting. Taking sex differences into account revealed that whereas stress over all increased wanting by 4.9 (CI=[0.56, 9.22], BF=69.2) points, stress increased drug wanting in only males. In males, wanting was 12.9 (CI=[6.72, 18.66], BF=Inf) percentage points higher in the stress compared to the control condition, and in females the differences was -3.1 (CI=[-9.1, 3.33], BF=0.2). The male-female differences of 16 (CI=[7.02, 24.49], BF=Inf) percentage points indicates that the effect of stress on drug wanting was higher in males than in females.

To investigate potential explanations for the sex difference in drug wanting we explored sex differences in the effect of drug on positive and negative affect, embarrassment/selfconsciousness and drug-side. This exploratory analysis showed that while females reported greater side effects (no differences for other affect variables), including drug effects on these potential mediators did not abolish the sex differences. This suggests that sex differences of drug effects on affect and side effects cannot fully explain the sex differences in drug wanting.

5.2 Effect of drug on post-stress induction affect

The following plots show changes in affect from pre drug administration (stage 4) to post drug administration (stage 5). The effect of drug is encoded in different pre-post slopes for the placebo and oxycodone conditions and can very by sex.

my_data.QPS = 
  my_data.Q %>% 
  .[, Stage := factor(stage.names[stage], levels = stage.names)] %>% 
  .[stage %in% c(4,5) & 
              item_name %in% c(negaffect_items,
                               posaffect_items,
                               embselfcn_items)] %>% 
  .[, affect := ifelse(item_name %in% negaffect_items,
                       "negative",
                       ifelse(item_name %in% posaffect_items,
                              "positive","emb_selfc"))] %>% 
  coarsen(var = "response",
          bins = 20,
          range = c(0,100))

my_data.QPS_stats = 
  my_data.QPS %>% 
  .[, .(response = mean(response),
        se = sd(response)/sqrt(.N)),
    by = .(Drug, stage, affect, Sex)] %>% 
  .[, `:=`(lower = response-2*se,
           upper = response+2*se)] 

my_data.QPS_stats %>% 
  ggplot(aes(x = stage, y = response, color = Drug,
             ymin = lower, ymax = upper, shape = Sex, lty = Sex)) +
  geom_point(position = pd.2) + 
  geom_line(aes(group = factor(Drug):factor(Sex)), position = pd.2) +
  geom_linerange(position = pd.2) + 
  facet_wrap(~affect) + 
  scale_x_continuous(breaks = c(3.75,4.75), labels = c("state\ninduction","dose1"), limits = c(3.5,5.25)) +
  scale_col_drug

The means already show clearly affect did not develop more positively in the oxycodone condition compared to the placebo condition.

We can make the same comparison for side effects:

my_data.QSE = 
  my_data.Q %>% 
  .[stage %in% c(4,5) & 
              item_name %in% sideffect_items] %>% 
  coarsen(var = "response",
          bins = 20,
          range = c(0,100))

my_data.QSE_stats = 
  my_data.QSE %>% 
  .[, .(response = mean(response),
        se = sd(response)/sqrt(.N)),
    by = .(Drug, stage, item_name, Sex)] %>% 
  .[, `:=`(lower = response-2*se,
           upper = response+2*se)] 

my_data.QSE_stats %>% 
    ggplot(aes(x = stage, y = response, color = Drug,
               ymin = lower, ymax = upper, shape = Sex, lty = Sex)) +
    geom_point(position = pd.2) + 
    geom_line(aes(group = factor(Drug):factor(Sex)), position = pd.2) +
    geom_linerange(position = pd.2) + 
    facet_wrap(~item_name) + 
  scale_x_continuous(breaks = c(3.75,4.75), labels = c("state\ninduction","dose1"), limits = c(3.5,5.25)) +
  scale_col_drug

This figure suggests that there was a greater increase in negative drug-related emotions in women than in men.

5.3 OLR with sex effects for placebo condition

my_data_beh_placebo = 
  copy(my_data_beh.all) %>% 
  .[Drug == "placebo"] %>% 
  .[effect_obtained > 125, effect_obtained := 125] %>% 
  coarsen(
    var = "effect_obtained",
    bins = n_cat,
    range = c(0,125))
analysis_formula = 
  formula(effect_obtained.c ~ state.condition:Sex + Sex + (1 + state.condition | participant) + (1|session))

We fit the model for effects obtained

options(mc.cores = 4)
fn = "brmsP1/obtained/placebo.Rdata"
if (file.exists(fn)) {
  load(fn)
} else {
  my_prior = 
    get_acat_thresh_priors(
      my_data_beh_placebo$effect_obtained.c,
      make_prior_b(prior_sd) + 
        prior(normal(0,2), class = "sd"))
  fit = brm(
    analysis_formula,
    family = acat(),
    data = my_data_beh_placebo,
    prior = my_prior,
    backend = "cmdstanr",
    iter = 2000)
  save(fit,file = fn)
}

And we make a figure as we did for the oxycodone condition:

plot.stats_table_obtained_placebo = 
  plot_results_origscale(
    fit,
    my_data_beh_placebo,
    outcome.var = "effect_obtained.c")

ggsave(plot.stats_table_obtained_placebo$plot,
       filename = "figures/EOplac.png", 
       width = 15,height = 15, units = "cm", device = cairo_pdf)

5.4 Stated desired effect

Prepare data

my_data_beh = 
  merge(my_data_beh,
        cheaters,
        by = c("participant","Drug","state.condition","session"),
        all.x = TRUE)
my_data_beh = 
  coarsen(my_data_beh,
          var = "target",
          bins = n_cat,
          range = c(0,125))

my_data_beh[is.na(cheater)] %>% 
  ggplot(aes(x = target, fill = state.condition)) + 
  geom_histogram(alpha = .5, position = "identity") + 
  scale_col_stress + 
  facet_grid(Sex~Drug) 

analysis_formula = 
  formula(target.c ~ cs(state.condition*Sex) + (1 + state.condition | participant) + (1|session))
fn = "brmsP1/desired.Rdata"
if (file.exists(fn)) {
  load(fn)
} else {
  my_prior = 
      prior(normal(0,2), class = b) + 
      prior(normal(0,1), class = Intercept)
  fit = brm(
    analysis_formula,
    family = acat(),
    data = my_data_beh[is.na(cheater)],
    prior = my_prior,
    iter = 5000,
    backend = "cmdstanr")
  save(fit,file = fn)
}

Make the plots

plot.stats_table_target = 
  plot_results_origscale(
    fit,
    my_data_beh[is.na(cheater)],
    outcome.var = "target.c")

p1 = plot.stats_table_target$p1 + ylab("% desired effect (0-125)")
yrange = layer_scales(p1)$y$range$range
p1 = p1 + scale_y_continuous(limits = yrange) +
  guides(y = guide_axis_truncated(
    trunc_lower = yrange[1]+diff(yrange)/500)) +
  coord_cartesian(ylim = c(67,108.5)) +
  xlab("") +
  theme(plot.title = element_text(hjust = -.05))

p2 = plot.stats_table_target$p2

px = 
  (p1 | p2) + 
  plot_layout(widths = c(4,1),) +
  plot_annotation(tag_levels = 'A')
ggsave(px,
        filename = "figures/Effect_Desired.png", 
        width = 20,height = 15, units = "cm")
ggsave(px,
        filename = "figures/Effect_Desired.pdf", 
        width = 20,height = 15, units = "cm")

Here are the estimated responses:

my_caption = "Model-estimated stated desired effect by sex and state condition."
my_footnote = "Numbers are means and CIs: Upper and lower bound of 95% credible intervals"
plot.stats_table_target$simple.stats %>% 
  kable(digits = 2, caption = my_caption) %>% 
  kable_styling() %>% 
  add_footnote(my_footnote, notation="none")

Table 5.9: Model-estimated stated desired effect by sex and state condition.

Sex	control	stress
men	90 (85, 95)	100 (95, 104)
total	84 (80, 88)	88 (85, 91)
women	78 (74, 83)	77 (72, 81)
Numbers are means and CIs: Upper and lower bound of 95% credible intervals

Here are the estimated effects:

“Effect of stress on desired effect in the self administration task”

5.5 Take drug again

Prepare data

my_data.Q.take_again = 
  my_data.Q %>% 
  .[Drug == "oxycodone" & item_name == "drug_again"] %>% 
  .[, .(participant,session,Sex,Drug,state.condition,stage, response)] %>% 
  coarsen(var = "response",
          bins = 20,
          range = c(0,100))

my_data.Q %>% 
  .[item_name == "drug_again"] %>% 
  ggplot(aes(x = response, fill = state.condition)) + 
  geom_histogram(alpha = .5, position = "identity") + 
  scale_col_stress + 
  facet_grid(Sex~Drug)

Figure 5.10: Willingess to take the drug again by gender and state condition.

Estimate the model (only with data from the oxycodon condition).

analysis_formula = 
  formula(response.c ~ cs(state.condition*Sex) + (1 + state.condition | participant) + (1|session))
options(mc.cores = 4)

fn = "brmsP1/again.Rdata"
if (file.exists(fn)) {
  load(fn)
} else {
  fit = brm(
    analysis_formula,
    family = acat(),
    data = my_data.Q.take_again,
    prior = my_prior,
    iter = 5000,
    backend = "cmdstanr")
  save(fit,file = fn)
}

Make the plots

plot.stats_table_again = 
  plot_results_origscale(
    fit,
    my_data.Q.take_again,
    outcome.var = "response.c")

Here are the estimated responses:

my_caption = "Model-estimated 'take again' responses by sex and state condition."
my_footnote = "Numbers are means and CIs: Upper and lower bound of 95% credible intervals"
plot.stats_table_again$simple.stats %>% 
  kable(digits = 2, caption = my_caption) %>% 
  kable_styling() %>% 
  add_footnote(my_footnote, notation="none")

Table 5.10: Model-estimated ‘take again’ responses by sex and state condition.

Sex	control	stress
men	31 (29, 33)	30 (27, 33)
total	28 (27, 30)	29 (27, 32)
women	26 (23, 29)	28 (26, 30)
Numbers are means and CIs: Upper and lower bound of 95% credible intervals

Here are the estimated effects:

cpt = "Effects of stress on 'take drug again'"
plot.stats_table_again$stats %>% 
  .[, CI := paste0("(",round(lower),", ", round(upper),")")] %>% 
  .[, c(1,3,8,7,6)] %>% 
  .[, `P > 0` := round(`P > 0`,2)] %>%
  .[, mean := round(mean)] %>%
  kable(caption = cpt) %>% 
  kable_styling(full_width = FALSE)

Table 5.11: Effects of stress on ‘take drug again’

Sex	mean	CI	BF	P > 0
men	-1	(-5, 2)	0.6	0.38
women	2	(-2, 5)	5.7	0.85
total	1	(-2, 3)	2.0	0.67
sex diff: m-f	-3	(-8, 2)	0.2	0.18

5.6 Plot for stated dersired effect and take drug again

Put plots together:

p = 
  (plot.stats_table_obtained$p1 + xlab("") + ylim(0,120) +
     scale_alpha(guide = 'none') +
     guides(y = guide_axis_truncated(trunc_lower = 0, trunc_upper = 120)) +
     theme(legend.position = "none", text = element_text(size = 10)) + 
     ylab("Drug self administration obtained (% of 1st)")) |
  (plot.stats_table_target$p1 + xlab("") + ylim(0,120) +
     guides(y = guide_axis_truncated(trunc_lower = 0, trunc_upper = 120)) +
     scale_alpha(guide = 'none') +
     theme(legend.position = "bottom", text = element_text(size = 10)) + 
     ylab("Drug self administration target (% of 1st)")) | 
  (plot.stats_table_again$p1 + ylim(0,100) +
     guides(y = guide_axis_truncated(trunc_lower = 0, trunc_upper = 100)) +
     scale_alpha(guide = 'none') +
     theme(legend.position = "none", text = element_text(size = 10)) + 
     ylab("Drug take again (%)")) +
  plot_layout(guides = "collect") +
  plot_annotation(tag_levels = 'A')
ggsave(p,filename = "figures/MR.png",width = 20,height = 10, units = "cm")

Figure X: Effect of stress on desired effect, obtained effect and rating to “take the drug again”.

6 Distribution of the number of button presses:

my_data_beh.all[, session := factor(session)]

my_data_beh.all = 
  merge(my_data_beh.all,
        cheaters,
        by = c("participant","Drug","state.condition","session"),
        all.x = TRUE) %>% 
  .[, n.presses.s := n.presses-mean(n.presses), by = .(participant)]


my_data_beh.all[, `target effect`  := ifelse(target > 50, "target > 50%","target < 50%")]

hist.presses = 
  my_data_beh.all[is.na(cheater)] %>% 
  ggplot(aes(x = n.presses, fill = Drug)) + 
  geom_histogram(position = "identity", alpha = .5, bins = 15) + 
  scale_fill_manual(values = palette()[c(2,4)]) +
  facet_grid(Sex~`target effect`) + 
  xlab("Number of presses") + 
  ylab("Number of participants") + 
  scale_col_drug

hist.time2target = 
  my_data_beh.all[is.na(cheater)] %>% 
  ggplot(aes(x = first.time.target, fill = Drug)) + 
  geom_histogram(position = "identity", alpha = .5, bins = 15) + 
  scale_fill_manual(values = palette()[c(2,4)]) +
   facet_grid(Sex~`target effect`) + 
  xlab("Time to target") + 
  ylab("Number of participants") + 
  scale_col_drug

#hist.presses / hist.time2target

my_data_beh.all %>% 
  .[, .(m = round(mean(n.presses)), sd = round(sd(n.presses))),
   by = .(Drug,state.condition,Sex)] %>% 
  kable(caption = "Number of button presses by druig and stress condition and sex") %>% 
  kable_styling(full_width = FALSE)

Table 6.1: Number of button presses by druig and stress condition and sex

Drug	state.condition	Sex	m	sd
oxycodone	control	women	613	101
oxycodone	stress	women	622	64
placebo	control	women	579	178
placebo	stress	women	615	138
oxycodone	control	men	639	161
oxycodone	stress	men	669	163
placebo	control	men	656	160
placebo	stress	men	670	135

# my_data_beh.all %>% 
#   ggplot(aes(x = state.condition, y = n.presses, color = Drug)) + 
#   stat_pointinterval(position = position_dodge(.2),  .width = c(0.5, 0.9)) + 
#   scale_color_manual(values = palette()[c(2,4)]) +
#   facet_wrap(~Sex) + 
#   scale_col_drug

7 Stress relief

To test if oxycodone buffers the effect of stress, we compare if the stress ratings increase more from pre to post state induction in the placebo condition compared to the oxycodone condition.

Here we plot the available data points:

my_data.Q.stressed %>% plot_stages()

The last pre state induction and pre drug time point is “BL 3” and the first post state induction and post drug time point is “dose 1”. Hence, we compare if the stress ratings increase more from “BL 3” to “dose 1” in the placebo condition compared to the oxycodone condition.

my_data.Q.stressdrug = select_data("stressed")
my_stages = c(3,5)
the_data = make_the_data(my_data.Q.stressdrug, "stressdrug")[, outcome := "stress"]

p1 = the_data %>% 
  .[, .(response = mean(response)), 
    by = .(participant, session, Drug, state.condition, stage)] %>% 
  ggplot(aes(x = stage, y = response, color = Drug,
             group = factor(participant):factor(session))) + 
  geom_point(alpha = .5) + 
  facet_wrap(~state.condition) + 
  geom_line(alpha = .5) + 
  scale_col_drug 

pdata = 
  the_data %>% 
  # average the pre stages BL1, BL2, BL3
  .[, .(response = mean(response)), 
    by = .(participant, Drug, state.condition, stage, Sex)] %>% 
  # calculate state induction effect as post-pre difference
  dcast(participant  + Drug + state.condition + Sex ~ stage, value.var = "response") %>%
  .[, induction_effect := post.stressdrug - pre.stressdrug] %>% 
  # calculate stress effect as difference of post-pre difference in stress & control conditions
  dcast(participant + Drug + Sex ~ state.condition, value.var = "induction_effect") %>% 
  .[, stress_effect := stress - control]

pdata.stats = 
  seWithin(pdata, "stress_effect",
           withinvars = c("Drug"),
           betweenvars = c("Sex"),
           idvar = "participant") %>% 
  data.table() %>% 
  .[, `:=`(lower=stress_effect-ci,
           upper=stress_effect+ci)]

p2 = pdata.stats %>% 
  ggplot(aes(x = Drug, y = stress_effect, fill = Drug)) + 
  geom_bar(stat = "identity", alpha = .5) + 
  geom_errorbar(aes(ymin = lower, ymax = upper), width = .25, size = 1) + 
  facet_wrap(~Sex) + 
  geom_beeswarm(data = pdata, aes(color = Drug)) +
  scale_col_drug +
  theme(legend.position = "none")

p3 = pdata  %>% 
  # calculate drug x stress interaction effect as difference of stress effects between drug conditions
  dcast(participant + Sex ~ Drug, value.var = "stress_effect") %>% 
  .[, stress_x_drug := oxycodone - placebo] %>% 
  ggplot(aes(x = stress_x_drug, fill=Sex)) + 
  geom_histogram(alpha = .5, position = "identity", bins = 20) +
  scale_col_sex +
  theme(legend.position = "right") 

p1 / (p2 | p3 ) + plot_annotation(tag_levels = 'A')

Figure 7.1: A: Stress ratings from pre to post state induction & drug administration. B: Effect of stress induction in the raw data measured after drug adminstration. Stress effects are calculated as the change in stress ratings from pre to post stress induction and drug administration. C: Distribution of the interaction of stress effect (pre-post induction) with drug condition (directly calculated from the data).

We coarsen the outcome variable in preparation of the ordinal regression model

the_data = the_data %>% 
  coarsen(var = "response", bins = 20, range = c(0,100))

stress_drug_formula = stress_drug_formula = response.c ~ induced.state:Sex:Drug + (1|session) + (1|participant/stage) + (induced.state:Sex:Drug|participant)

stress_drug_formula = stress_drug_formula = response.c ~ induced.state:Sex:Drug + (1|session) + (1|participant/stage) + (0+induced.state:Sex:Drug|participant)

We fit the model

fn = "brmsP1/stress_relief_acat.Rdata"
if (file.exists(fn)) {
  load(fn)
} else {
  fit = fit_ordinal_model(the_data, stress_drug_formula, fn, family = acat())
}
#pp_check(fit, ndraws = 2000, type = "bars")

Click to show sample diagnostics

Table 7.1: Summary of model paramters

parameter	mean	sd	q5	q95	rhat	ess_bulk	ess_tail
Threshold[1]	1.0	0.73	-0.2	2.2	1.002	784	1312
Threshold[2]	0.5	0.73	-0.7	1.7	1.002	803	1438
Threshold[3]	0.5	0.74	-0.8	1.7	1.002	821	1259
Threshold[4]	0.9	0.75	-0.3	2.2	1.002	839	1484
Threshold[5]	0.5	0.76	-0.8	1.7	1.002	869	1668
Threshold[6]	0.6	0.78	-0.7	1.9	1.002	872	1516
Threshold[7]	0.6	0.79	-0.6	1.9	1.003	934	1559
Threshold[8]	0.9	0.82	-0.4	2.2	1.002	910	1587
Threshold[9]	0.2	0.82	-1.1	1.5	1.002	976	1753
Threshold[10]	1.4	0.86	-0.1	2.8	1.002	1017	1901
Threshold[11]	0.9	0.92	-0.6	2.4	1.002	1146	2017
Threshold[12]	1.2	1.03	-0.5	2.9	1.002	1273	2217
Threshold[13]	0.1	1.00	-1.6	1.7	1.001	1368	2118
Threshold[14]	2.5	1.20	0.6	4.6	1.001	1462	2341
Threshold[15]	0.6	1.39	-1.6	3.0	1.001	2255	2693
b_induced.statecontrol:Sexman:Drugoxycodone	-0.1	0.70	-1.3	1.0	1.002	811	1459
b_induced.statestress:Sexman:Drugoxycodone	0.0	0.72	-1.2	1.1	1.002	808	1361
b_induced.statecontrol:Sexwoman:Drugoxycodone	0.0	0.70	-1.2	1.1	1.003	802	1433
b_induced.statestress:Sexwoman:Drugoxycodone	0.2	0.71	-1.0	1.3	1.003	809	1428
b_induced.statecontrol:Sexman:Drugplacebo	-0.1	0.71	-1.3	1.0	1.003	799	1454
b_induced.statestress:Sexman:Drugplacebo	0.0	0.71	-1.1	1.2	1.003	814	1438
b_induced.statecontrol:Sexwoman:Drugplacebo	-0.1	0.71	-1.2	1.1	1.003	804	1415
b_induced.statestress:Sexwoman:Drugplacebo	0.2	0.71	-0.9	1.4	1.004	791	1459
sd_participant__Intercept	0.3	0.07	0.2	0.4	1.006	826	799
sd_participant__induced.statecontrol:Sexman:Drugoxycodone	0.3	0.15	0.1	0.6	1.001	711	722
sd_participant__induced.statestress:Sexman:Drugoxycodone	0.2	0.18	0.0	0.6	1.003	1489	1962
sd_participant__induced.statecontrol:Sexwoman:Drugoxycodone	0.3	0.12	0.1	0.5	1.002	957	781
sd_participant__induced.statestress:Sexwoman:Drugoxycodone	0.3	0.17	0.0	0.6	1.006	784	1154
sd_participant__induced.statecontrol:Sexman:Drugplacebo	0.2	0.14	0.0	0.5	1.007	578	1194
sd_participant__induced.statestress:Sexman:Drugplacebo	0.2	0.14	0.0	0.5	1.007	1629	2090
sd_participant__induced.statecontrol:Sexwoman:Drugplacebo	0.3	0.14	0.1	0.6	1.003	686	642
sd_participant__induced.statestress:Sexwoman:Drugplacebo	0.2	0.16	0.0	0.5	1.002	909	1888
sd_participant:stage__Intercept	0.1	0.06	0.0	0.2	1.004	621	1148
sd_session__Intercept	0.3	0.21	0.1	0.6	1.008	1041	1451
model: response.c ~ induced.state:Sex:Drug + (1 | session) + (1 | participant/stage) + (0 + induced.state:Sex:Drug | participant) b_ = fix effect parameters sd_ standard deviation of hierarchical parameters ess = effective sample size

We calculate the same contrasts as above:

• stress effect under oxycodone: (((pre-post|stress,oxy)-(pre-post|control,oxy)) minus
• stress effect under placebo: ((pre-post|stress,plac)-(pre-post|control,plac))`)

this time on from the posterior predictions:

fn = "brmsP1/stress_relief_contrasts.Rdata"
if (!file.exists(fn)) {
  postpre = unique(the_data$stage) %>% as.character() %>% sort()
  new_data = do.call(rbind,lapply(1:4, function(x) data.table(fit$data)[, session := x]))
  epred = posterior_epred(fit, newdata = new_data)
  obs.ratings = the_data[,.(r = mean(response)), by = .(response.c)][,r] %>% sort()
  epred.r = 
    apply(epred, 1, function(a) apply(a, 1, function(p) sum(p*obs.ratings)))
  
  epred.data = 
    new_data %>% 
    .[, .(induced.state, stage, Sex, participant, Drug)] %>% 
    .[, state.condition := ifelse(any(induced.state == "stress"),"stress","control"), 
      by = .(participant,Drug,induced.state)] %>% 
    .[, induced.state := NULL]
  
  pp = 
    cbind(epred.data, epred.r) %>% 
    melt(id.vars = names(epred.data), variable.name = "iter") %>% 
    .[, .(value = mean(value)), by = .(Sex,Drug,state.condition, iter)] %>% 
    .[, Sex := gsub("man","men",Sex)]
  
  
  # posterior distribution of contrasts, i.e.
  # main effects in men and women + their average and difference
  # main effects in men and women
  px = pp %>% 
    dcast(iter + Sex ~ Drug + state.condition) %>% 
    .[, stress_x_drug := 
        (oxycodone_stress - oxycodone_control) - 
        (placebo_stress - placebo_control)] %>% 
    .[, stress := 
        ((oxycodone_stress - oxycodone_control) + 
        (placebo_stress - placebo_control))/2] %>% 
    .[, .(iter, Sex, stress_x_drug, stress)]
  
  px = rbind(
    px,
    # average
    px[, .(stress_x_drug = mean(stress_x_drug),
           stress = mean(stress)),
       by = .(iter)][, Sex := "total"],
    # sex difference
    px %>% dcast(iter ~ Sex, value.var = c("stress","stress_x_drug")) %>% 
      .[, `:=`(stress_x_drug = stress_x_drug_men-stress_x_drug_women,
               stress = stress_men-stress_women)] %>% 
      .[, Sex := "Sex difference"] %>% 
      .[, .(iter, stress, stress_x_drug, Sex)]
  ) %>% 
    .[, Sex := factor(Sex, levels = c("men","women","Sex difference","total"))]
  save(px, pp, file = fn)
} else {
  load(fn)
}

Consistent with the simple stress effects seen in the raw data, we a credible stress effect is visible only in women.

p1 = pp %>% 
  dcast(iter + Sex + Drug ~ state.condition) %>% 
  .[, stress_effect := stress-control] %>% 
  ggplot(aes(x = Sex, y = stress_effect, color = Drug)) +
  stat_pointinterval(position = position_dodge(.2)) + 
  geom_hline(yintercept = 0, col = "red", lty = 2) + 
  scale_col_drug 

p2 = px %>% 
  ggplot(aes(x = stress_x_drug, fill = Sex)) + 
  geom_vline(xintercept = 0) + 
  geom_density(alpha = .5, color = NA) + 
  facet_wrap(~Sex) + 
  scale_col_sex.b

(p1 | p2) + 
  plot_layout(widths = c(2,3)) + 
  plot_annotation(tag_levels = "A")

Figure 7.2: A: Model-predicted effects of stress induction on ‘stressed’ ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Stress by Drug interaction in women, men, and sex difference (women-men)

We found no reliable effect of drug on stress. Oxycodone reduced stress induction induced increase of stress rating from pre state induction to post drug administration on average by NaN (CI=[NA, NA]; BFNA; P(O>P)NA). The interaction effect was -1 (CI=[-7, 6]; BF=0.7; P(O>P)=0.4) in women and -1 (CI=[-6, 5]; BF=0.6; P(O>P)=0.38) in men.

8 Buffering effects for reminders

Prepare the data:

my_data.Q_buffer = 
  my_data.Q %>% 
  .[stage %in% c(5,6) & 
      item_name %in% c(negaffect_items,
                       posaffect_items,
                       embselfcn_items,
                       "stressed")] %>% 
  .[, stress_reminder := ifelse(stage %in% c(5,7) & state.condition == "stress","pre","post")] %>% 
  .[, stress_reminder := factor(stress_reminder, levels = c("pre","post"))] %>% 
  .[, reminder := ifelse(stage %in% c(5,7),"pre","post")] %>% 
  .[, reminder := factor(reminder, levels = c("pre","post"))] %>% 
  .[, affect := ifelse(item_name %in% negaffect_items,
                       "negative",
                       ifelse(item_name %in% posaffect_items,
                              "positive","emb_selfc"))] %>% 
  .[item_name == "stressed", affect := "stressed"] %>% 
  coarsen(var = "response",
          bins = 20,
          range = c(0,100))

Here is a plot of the (nearly) raw data.

my_data.Q_buffer_stats = 
  my_data.Q_buffer %>% 
  .[, .(response = mean(response),
        se = sd(response)/sqrt(.N)),
    by = .(Drug, reminder, item_name, Sex, state.condition)] %>% 
  .[, `:=`(lower = response-2*se,
           upper = response+2*se)] %>% 
  .[, grp := paste(Drug,state.condition)]

my_data.Q_buffer_stats %>% 
  ggplot(aes(x = reminder, y = response, color = Drug, lty = state.condition, 
             shape = state.condition, group = grp,
             ymin = lower, ymax = upper)) +
  geom_point(position = pd.2) + 
  geom_line(aes(group = grp), position = pd.2) +
  geom_linerange(position = pd.2) + 
  facet_grid(item_name~Sex, scales = "free_y") +
  scale_col_drug

Figure 8.1: Mean ratings for stress, embarrassment and selfconciousness pre and post reminders by drug, state condition and sex.

pdata = 
  my_data.Q_buffer %>% 
  # average the pre stages BL1, BL2, BL3
  .[, .(response = mean(response)), 
    by = .(participant, Drug, state.condition, stage, Sex,item_name )] %>% 
  # calculate state induction effect as post-pre difference
  dcast(participant  + Drug + state.condition + Sex + item_name ~ stage, value.var = "response") %>%
  .[, induction_effect := `6` - `5`] %>% 
  # calculate stress effect as difference of post-pre difference in stress & control conditions
  dcast(participant + Drug + Sex + item_name ~ state.condition, value.var = "induction_effect") %>% 
  .[, stress_effect := stress - control]

pdata.stats = 
  seWithin(pdata, "stress_effect",
           withinvars = c("Drug","item_name"),
           betweenvars = c("Sex"),
           idvar = "participant") %>% 
  data.table() %>% 
  .[, `:=`(lower=stress_effect-ci,
           upper=stress_effect+ci)]

pdata.stats %>% 
  ggplot(aes(x = Drug, y = stress_effect, fill = Drug)) + 
  geom_bar(stat = "identity", alpha = .5) + 
  geom_errorbar(aes(ymin = lower, ymax = upper), width = .25, linewidth = 1) + 
  facet_grid(item_name~Sex) + 
  geom_quasirandom(data = pdata, aes(color = Drug)) +
  scale_col_drug +
  ylim(-50,50)

Figure 8.2: Effect of stress reminders in the raw data measured after drug adminstration. Reminder effects are calculated as the change in stress ratings from pre to post stress induction and drug administration.

pdata %>% 
  .[, stress_effect := stress - control] %>% 
  # calculate drug x stress interaction effect as difference of stress effects between drug conditions
  dcast(participant + item_name + Sex ~ Drug, value.var = "stress_effect") %>% 
  .[, stress_x_drug := oxycodone - placebo] %>% 
  ggplot(aes(x = stress_x_drug, fill=Sex)) + 
  geom_histogram(alpha = .5, position = "identity", bins = 20) + 
  scale_col_sex + 
  facet_wrap(~item_name)

Figure 8.3: Distribution of the interaction of stress effect (pre-post reminder) with drug condition (directly calculated from the data).

8.1 Analysis approach

The analysis approach here is that we get posterior distributions for the ratings in the 8 cells obtained by crossing:

• stress vs control
• pre vs post reminder
• drug vs placebo
• males vs females

Then we can calculate the desired contrasts by first calculating the pre to post change for 4 conditions from stress x drug. - pre_post_placebo_control = Drugplacebo:Stresscontrol:reminderpre - Drugplacebo:Stresscontrol:reminderpost - pre_post_placebo_stress = Drugplacebo:Stressstress:reminderpre - Drugplacebo:Stressstress:reminderpost - pre_post_placebo_control = Drugplacebo:Stresscontrol:reminderpre - Drugplacebo:Stresscontrol:reminderpost - pre_post_placebo_stress = Drugplacebo:Stressstress:reminderpre - Drugplacebo:Stressstress:reminderpost

With those, we can calculate the stress reminder effects in the drug and placebo condition

• stress_effect_placebo = pre_post_placebo_stress - pre_post_placebo_control
• stress_effect_placebo = pre_post_placebo_stress - pre_post_placebo_control

And finally we can calculate if the stress-reminder effect is smaller in the the placebo compared to the placebo condition

• stress_x_drug = stress_effect_placebo - stress_effect_placebo

This description does not include sex effects, but these effects can be calculates separately by sex.

Note that given this specification of the contrasts, a buffer effect where post-reminder ratings of negative items like stress are lower under oxycodone compared to placebo will result in negative values for the stress_x_drug interaction effect.

Here is the analysis model:

analysis_formula = 
  response.c ~  Sex:Drug:Stress:reminder + (1|session) + (1 + Drug:Stress:reminder | participant) 

The analysis for manipulation checks showed that that the stress effects of the first reminder were clearer, especially for feeling stressed. Therefore, the analysis of buffer effects focuses on the first reminder.

8.2 Buffering of stress

We focus on buffering of stress for the first reminder.

f = "brmsP1/buffer/stressed.Rdata"
the_data = my_data.Q_buffer[item_name == "stressed"]
if (file.exists(f)) {
  load(f)
} else {
  fit = fit_ordinal_model(
    the_data, analysis_formula, f, family = acat())
}

Click to show sample diagnostics

Table 8.1: Summary of model paramters

parameter	mean	sd	q5	q95	rhat	ess_bulk	ess_tail
Threshold[1]	-2.0	0.97	-3.6	-0.4	1.014	394	409
Threshold[2]	-1.7	0.88	-3.2	-0.3	1.013	452	521
Threshold[3]	-0.8	0.84	-2.2	0.6	1.013	500	755
Threshold[4]	-0.6	0.84	-2.0	0.8	1.012	566	801
Threshold[5]	-0.3	0.84	-1.7	1.1	1.011	624	938
Threshold[6]	0.0	0.84	-1.4	1.4	1.012	646	1000
Threshold[7]	0.5	0.86	-0.9	1.9	1.010	770	1494
Threshold[8]	0.3	0.87	-1.2	1.7	1.010	726	1568
Threshold[9]	0.8	0.86	-0.7	2.2	1.008	828	1672
Threshold[10]	1.6	0.93	0.0	3.1	1.008	987	1638
Threshold[11]	1.4	0.99	-0.2	3.1	1.006	1091	2064
Threshold[12]	3.0	1.22	1.0	5.0	1.003	1511	2240
Threshold[13]	-0.5	1.17	-2.4	1.4	1.006	1385	2313
Threshold[14]	3.4	1.08	1.6	5.2	1.006	1327	2325
Threshold[15]	1.9	1.16	-0.1	3.8	1.004	1552	2328
Threshold[16]	3.1	1.17	1.2	5.1	1.004	1440	2207
Threshold[17]	3.5	1.27	1.4	5.6	1.004	1724	2671
b_Sexman:Drugoxycodone:Stresscontrol:reminderpre	-1.1	0.72	-2.3	0.1	1.007	894	1487
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpre	-0.5	0.68	-1.7	0.6	1.008	883	1976
b_Sexman:Drugplacebo:Stresscontrol:reminderpre	-0.8	0.63	-1.9	0.2	1.007	774	1486
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpre	-0.2	0.62	-1.2	0.8	1.007	869	1573
b_Sexman:Drugoxycodone:Stressstress:reminderpre	-0.2	0.66	-1.3	0.9	1.004	938	1380
b_Sexwoman:Drugoxycodone:Stressstress:reminderpre	0.7	0.65	-0.4	1.7	1.007	905	1500
b_Sexman:Drugplacebo:Stressstress:reminderpre	0.2	0.60	-0.8	1.2	1.004	895	1606
b_Sexwoman:Drugplacebo:Stressstress:reminderpre	0.8	0.60	-0.2	1.8	1.005	857	1584
b_Sexman:Drugoxycodone:Stresscontrol:reminderpost	-1.0	0.77	-2.3	0.3	1.005	1162	1845
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpost	-0.6	0.75	-1.8	0.6	1.005	1081	2023
b_Sexman:Drugplacebo:Stresscontrol:reminderpost	-1.2	0.67	-2.3	-0.1	1.007	834	1387
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpost	-0.6	0.67	-1.8	0.4	1.005	943	1786
b_Sexman:Drugoxycodone:Stressstress:reminderpost	0.8	0.66	-0.3	1.9	1.007	882	1465
b_Sexwoman:Drugoxycodone:Stressstress:reminderpost	1.6	0.67	0.5	2.7	1.007	848	1549
b_Sexman:Drugplacebo:Stressstress:reminderpost	0.1	0.65	-1.0	1.2	1.004	973	1655
b_Sexwoman:Drugplacebo:Stressstress:reminderpost	1.9	0.66	0.8	3.0	1.007	831	1539
sd_participant__Intercept	1.5	0.32	1.0	2.0	1.003	697	1364
sd_participant__Drugoxycodone:Stresscontrol:reminderpre	1.9	0.48	1.2	2.8	1.003	606	1339
sd_participant__Drugplacebo:Stresscontrol:reminderpre	0.7	0.42	0.1	1.5	1.005	435	1106
sd_participant__Drugoxycodone:Stressstress:reminderpre	1.6	0.43	1.0	2.4	1.006	591	1150
sd_participant__Drugplacebo:Stressstress:reminderpre	1.0	0.39	0.4	1.7	1.013	410	566
sd_participant__Drugoxycodone:Stresscontrol:reminderpost	2.2	0.50	1.5	3.1	1.010	527	877
sd_participant__Drugplacebo:Stresscontrol:reminderpost	1.1	0.47	0.4	2.0	1.013	389	440
sd_participant__Drugoxycodone:Stressstress:reminderpost	2.0	0.46	1.3	2.8	1.003	435	720
sd_participant__Drugplacebo:Stressstress:reminderpost	1.5	0.44	0.8	2.3	1.009	380	535
sd_session__Intercept	2.6	0.93	1.2	4.2	1.003	632	405
model: response.c ~ Sex:Drug:Stress:reminder + (1 | session) + (1 + Drug:Stress:reminder | participant) b_ = fix effect parameters sd_ standard deviation of hierarchical parameters ess = effective sample size

p = plot_buffer_effect(fit,the_data,"stressed")
p$ppc

Figure 8.4: Posterior predictive check of pre-post differences. Points are differences directly calculate from the data, intervals show differences calculate from model parameters.

The following figure shows stress-drug interactions (negative values indicate a buffering effect). The data support a hypothesis of no effect of oxycodone on stress response if the credible interval falls in a region of practical equivalence of +/- 5 points.

((p$stress.plot + 
    ylab("stress effect on 'stressed' ratings.")) | p$buffer.plot) + 
  plot_annotation(tag_levels = "A") + 
  plot_layout(widths = c(2,3))

Figure 8.5: A: Model-predicted effects of the first stress reminder on ‘stressed’ ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Posterior distribution of the interaction of stress effect (pre-post reminder) on ‘stressed’ ratings with drug condition. Dashed vertical lines enclode a region of practical equivalence.

p$buffer.stats %>% 
  kable(caption = "Buffering effect of oxycodone for 'stressed' ratings after reminders.") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic()

Table 8.2: Buffering effect of oxycodone for ‘stressed’ ratings after reminders.

Sex	stress	stress in placebo	stress in oxycodone	buffer
men	4 [2, 7]; >.999	4 [0, 7]; 0.9775	5 [1, 9]; 0.997	2 [-4, 7]; 0.72725
women	12 [9, 15]; >.999	16 [12, 20]; >.999	8 [3, 12]; >.999	-8 [-14, -2]; 0.0035
total	8 [6, 10]; >.999	10 [7, 13]; >.999	6 [4, 9]; >.999	-3 [-7, 1]; 0.05
Sex difference	7 [3, 11]; >.999	12 [6, 18]; >.999	2 [-3, 8]; 0.7955	-10 [-18, -2]; 0.0095

8.3 Buffering of embarrassment

f = "brmsP1/buffer/emb.Rdata"
the_data = my_data.Q_buffer[item_name == "embarassed"]
if (file.exists(f)) {
  load(f)
} else {
  fit = fit_ordinal_model(
    the_data, analysis_formula, f, family = acat())
}

Click to show sample diagnostics

Table 8.3: Summary of model paramters

parameter	mean	sd	q5	q95	rhat	ess_bulk	ess_tail
Threshold[1]	-2.0	0.87	-3.4	-0.5	1.003	920	1447
Threshold[2]	-1.5	0.80	-2.8	-0.2	1.002	990	1480
Threshold[3]	-1.2	0.78	-2.5	0.1	1.001	1067	1540
Threshold[4]	-0.5	0.79	-1.8	0.8	1.002	1113	1589
Threshold[5]	-0.2	0.81	-1.5	1.2	1.002	1237	1720
Threshold[6]	0.2	0.84	-1.2	1.6	1.001	1265	1890
Threshold[7]	0.5	0.90	-0.9	2.0	1.002	1628	2181
Threshold[8]	0.5	0.94	-1.0	2.0	1.001	1682	1907
Threshold[9]	-0.5	0.91	-2.1	1.0	1.001	1517	2289
Threshold[10]	2.0	0.95	0.4	3.5	1.000	1745	2642
Threshold[11]	0.8	1.06	-0.9	2.6	1.000	1966	2063
Threshold[12]	0.3	1.00	-1.3	2.0	1.001	1901	2965
Threshold[13]	1.3	0.90	-0.1	2.8	1.000	1521	2469
Threshold[14]	1.3	0.88	-0.2	2.7	1.001	1576	2248
Threshold[15]	4.3	1.20	2.4	6.3	1.002	2299	2466
Threshold[16]	1.4	1.28	-0.6	3.5	1.000	2606	2651
Threshold[17]	2.2	1.15	0.3	4.1	1.002	2114	2831
Threshold[18]	2.5	1.08	0.8	4.3	1.000	1859	2689
b_Sexman:Drugoxycodone:Stresscontrol:reminderpre	-1.1	0.73	-2.4	0.0	1.001	1552	2509
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpre	-0.4	0.70	-1.6	0.8	1.001	1624	2183
b_Sexman:Drugplacebo:Stresscontrol:reminderpre	-0.2	0.71	-1.4	1.0	1.001	1723	2268
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpre	-1.3	0.81	-2.7	0.0	1.002	1760	1997
b_Sexman:Drugoxycodone:Stressstress:reminderpre	-0.3	0.73	-1.5	0.9	1.002	1590	2467
b_Sexwoman:Drugoxycodone:Stressstress:reminderpre	0.3	0.75	-1.0	1.5	1.001	1533	1957
b_Sexman:Drugplacebo:Stressstress:reminderpre	-0.4	0.73	-1.6	0.8	1.001	1805	2477
b_Sexwoman:Drugplacebo:Stressstress:reminderpre	1.3	0.68	0.1	2.4	1.001	1545	2003
b_Sexman:Drugoxycodone:Stresscontrol:reminderpost	-0.7	0.71	-1.9	0.4	1.002	1521	2407
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpost	-1.6	0.78	-2.9	-0.3	1.003	1421	2283
b_Sexman:Drugplacebo:Stresscontrol:reminderpost	-1.0	0.67	-2.1	0.1	1.001	1572	2120
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpost	-1.6	0.75	-2.8	-0.4	1.001	1528	2385
b_Sexman:Drugoxycodone:Stressstress:reminderpost	0.9	0.71	-0.2	2.1	1.001	1640	2468
b_Sexwoman:Drugoxycodone:Stressstress:reminderpost	1.6	0.72	0.4	2.8	1.001	1555	2074
b_Sexman:Drugplacebo:Stressstress:reminderpost	0.9	0.67	-0.1	2.0	1.002	1582	2393
b_Sexwoman:Drugplacebo:Stressstress:reminderpost	3.5	0.75	2.3	4.8	1.002	1254	2026
sd_participant__Intercept	1.9	0.39	1.4	2.6	1.002	1303	1937
sd_participant__Drugoxycodone:Stresscontrol:reminderpre	1.3	0.48	0.5	2.1	1.008	1050	965
sd_participant__Drugplacebo:Stresscontrol:reminderpre	1.7	0.45	1.0	2.5	1.002	1397	2582
sd_participant__Drugoxycodone:Stressstress:reminderpre	2.1	0.50	1.3	2.9	1.002	1306	2050
sd_participant__Drugplacebo:Stressstress:reminderpre	1.8	0.45	1.1	2.5	1.010	783	1111
sd_participant__Drugoxycodone:Stresscontrol:reminderpost	1.0	0.50	0.2	1.8	1.014	662	1056
sd_participant__Drugplacebo:Stresscontrol:reminderpost	0.6	0.44	0.1	1.5	1.007	816	1699
sd_participant__Drugoxycodone:Stressstress:reminderpost	2.1	0.50	1.4	3.0	1.002	1102	1792
sd_participant__Drugplacebo:Stressstress:reminderpost	1.8	0.44	1.1	2.6	1.015	745	1178
sd_session__Intercept	3.0	0.93	1.7	4.7	1.003	1435	1728
model: response.c ~ Sex:Drug:Stress:reminder + (1 | session) + (1 + Drug:Stress:reminder | participant) b_ = fix effect parameters sd_ standard deviation of hierarchical parameters ess = effective sample size

p = plot_buffer_effect(fit,the_data,"embarassed")
p$ppc

Figure 8.6: Posterior predictive check of pre-post differences. Points are differences directly calculate from the data, intervals show differences calculate from model parameters.

((p$stress.plot + 
    ylab("stress effect on 'embarassed' ratings.")) | p$buffer.plot) + 
  plot_annotation(tag_levels = "A") + 
  plot_layout(widths = c(2,3))

Figure 8.7: A: Model-predicted effects of the first stress reminder on ‘embarassed’ ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Posterior distribution of the interaction of stress effect (pre-post reminder) on ‘embarassed’ ratings with drug condition. Dashed vertical lines enclode a region of practical equivalence.

p$buffer.stats %>% 
  kable(caption = "Buffering effect of oxycodone for 'embarassed' ratings after reminders.") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic()

Table 8.4: Buffering effect of oxycodone for ‘embarassed’ ratings after reminders.

Sex	stress	stress in placebo	stress in oxycodone	buffer
men	9 [7, 12]; >.999	11 [7, 14]; >.999	7 [4, 11]; >.999	-3 [-8, 2]; 0.09525
women	22 [20, 25]; >.999	30 [25, 34]; >.999	15 [11, 19]; >.999	-15 [-20, -9]; <.001
total	16 [14, 18]; >.999	20 [17, 23]; >.999	11 [9, 14]; >.999	-9 [-13, -5]; <.001
Sex difference	13 [10, 17]; >.999	19 [13, 25]; >.999	8 [3, 13]; 0.9985	-11 [-19, -3]; 0.00225

8.4 Buffering of selfconciousness

f = "brmsP1/buffer/selfcn.Rdata"
the_data = my_data.Q_buffer[item_name == "selfconscious"]
if (file.exists(f)) {
  load(f)
} else {
  fit = fit_ordinal_model(
    the_data, analysis_formula, f, family = acat())
}

Click to show sample diagnostics

Table 8.5: Summary of model paramters

parameter	mean	sd	q5	q95	rhat	ess_bulk	ess_tail
Threshold[1]	-2.3	0.78	-3.5	-1.0	1.012	500	1004
Threshold[2]	-1.6	0.80	-3.0	-0.3	1.008	672	1286
Threshold[3]	-3.2	0.80	-4.5	-1.9	1.007	703	1569
Threshold[4]	-1.8	0.74	-3.0	-0.6	1.011	648	1189
Threshold[5]	-1.8	0.73	-3.0	-0.6	1.009	665	1429
Threshold[6]	-2.2	0.70	-3.4	-1.1	1.010	609	1349
Threshold[7]	-1.3	0.66	-2.4	-0.3	1.010	603	1289
Threshold[8]	-0.3	0.71	-1.5	0.9	1.007	707	1567
Threshold[9]	-2.5	0.67	-3.7	-1.5	1.008	655	1381
Threshold[10]	0.1	0.61	-0.9	1.0	1.008	581	1214
Threshold[11]	-0.1	0.63	-1.2	0.9	1.009	585	1271
Threshold[12]	-0.1	0.61	-1.2	0.8	1.006	578	1219
Threshold[13]	0.8	0.63	-0.2	1.8	1.006	571	1200
Threshold[14]	0.8	0.64	-0.3	1.8	1.006	542	1124
Threshold[15]	1.8	0.66	0.7	2.9	1.006	555	1182
Threshold[16]	1.2	0.69	0.1	2.4	1.002	522	1044
Threshold[17]	2.3	0.70	1.2	3.5	1.004	491	964
Threshold[18]	3.5	0.76	2.2	4.7	1.002	502	1011
Threshold[19]	2.7	0.83	1.3	4.0	1.003	450	909
b_Sexman:Drugoxycodone:Stresscontrol:reminderpre	0.0	0.64	-1.0	1.1	1.009	453	867
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpre	-0.4	0.61	-1.4	0.6	1.006	523	987
b_Sexman:Drugplacebo:Stresscontrol:reminderpre	0.0	0.62	-1.0	1.0	1.008	449	918
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpre	-0.2	0.59	-1.1	0.8	1.007	505	1061
b_Sexman:Drugoxycodone:Stressstress:reminderpre	0.2	0.63	-0.9	1.2	1.006	477	1138
b_Sexwoman:Drugoxycodone:Stressstress:reminderpre	-0.5	0.60	-1.5	0.5	1.005	575	1181
b_Sexman:Drugplacebo:Stressstress:reminderpre	0.4	0.63	-0.6	1.4	1.006	481	836
b_Sexwoman:Drugplacebo:Stressstress:reminderpre	0.1	0.60	-0.8	1.1	1.005	510	969
b_Sexman:Drugoxycodone:Stresscontrol:reminderpost	0.4	0.62	-0.7	1.4	1.007	427	867
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpost	-0.5	0.60	-1.5	0.4	1.006	510	1034
b_Sexman:Drugplacebo:Stresscontrol:reminderpost	-0.2	0.65	-1.3	0.9	1.007	437	904
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpost	-0.8	0.63	-1.8	0.3	1.006	544	1040
b_Sexman:Drugoxycodone:Stressstress:reminderpost	0.3	0.61	-0.7	1.3	1.006	495	926
b_Sexwoman:Drugoxycodone:Stressstress:reminderpost	-0.2	0.58	-1.1	0.8	1.006	550	1158
b_Sexman:Drugplacebo:Stressstress:reminderpost	0.5	0.64	-0.6	1.5	1.006	512	1014
b_Sexwoman:Drugplacebo:Stressstress:reminderpost	1.0	0.63	0.0	2.0	1.004	559	954
sd_participant__Intercept	2.4	0.41	1.8	3.1	1.010	406	923
sd_participant__Drugoxycodone:Stresscontrol:reminderpre	1.4	0.36	0.9	2.0	1.015	296	341
sd_participant__Drugplacebo:Stresscontrol:reminderpre	0.7	0.38	0.1	1.4	1.027	154	287
sd_participant__Drugoxycodone:Stressstress:reminderpre	1.5	0.33	1.0	2.1	1.009	364	871
sd_participant__Drugplacebo:Stressstress:reminderpre	0.9	0.34	0.4	1.5	1.010	150	355
sd_participant__Drugoxycodone:Stresscontrol:reminderpost	1.1	0.33	0.6	1.7	1.014	260	235
sd_participant__Drugplacebo:Stresscontrol:reminderpost	1.3	0.35	0.7	1.9	1.007	269	944
sd_participant__Drugoxycodone:Stressstress:reminderpost	1.2	0.30	0.8	1.8	1.018	368	772
sd_participant__Drugplacebo:Stressstress:reminderpost	1.6	0.36	1.1	2.2	1.004	259	861
sd_session__Intercept	0.1	0.16	0.0	0.4	1.001	1390	1909
model: response.c ~ Sex:Drug:Stress:reminder + (1 | session) + (1 + Drug:Stress:reminder | participant) b_ = fix effect parameters sd_ standard deviation of hierarchical parameters ess = effective sample size

p = plot_buffer_effect(fit,the_data,"selfconscious")
p$ppc

Figure 8.8: Posterior predictive check of pre-post differences. Points are differences directly calculate from the data, intervals show differences calculate from model parameters.

((p$stress.plot + 
    ylab("stress effect on 'selfconscious' ratings.")) | p$buffer.plot) + 
  plot_annotation(tag_levels = "A") + 
  plot_layout(widths = c(2,3))

Figure 8.9: A: Model-predicted effects of the first stress reminder on ‘selfconscious’ ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Posterior distribution of the interaction of stress effect (pre-post reminder) on ‘selfconscious’ ratings with drug condition. Dashed vertical lines enclode a region of practical equivalence.

p$buffer.stats %>% 
  kable(caption = "Buffering effect of oxycodone for 'selfconscious' ratings after reminders.") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic()

Table 8.6: Buffering effect of oxycodone for ‘selfconscious’ ratings after reminders.

Sex	stress	stress in placebo	stress in oxycodone	buffer
men	1 [-3, 5]; 0.73425	4 [-1, 10]; 0.94775	-2 [-7, 3]; 0.244	-6 [-13, 1]; 0.051
women	12 [8, 16]; >.999	18 [12, 24]; >.999	6 [0, 12]; 0.97675	-12 [-20, -4]; 0.0025
total	7 [4, 9]; >.999	11 [7, 15]; >.999	2 [-2, 6]; 0.85225	-9 [-14, -4]; <.001
Sex difference	11 [5, 17]; >.999	14 [6, 22]; >.999	8 [-0, 16]; 0.97175	-6 [-17, 5]; 0.14725

8.5 Buffering of stress, embarrassement & selfconciousness

Given that the data patterns for these three items are expected to be similar, we can obtain more stable estimates of a buffer effect by modeling all items jointly.

f = "brmsP1/buffer/embselfstress.Rdata"
the_data = my_data.Q_buffer[item_name %in% c(embselfcn_items,"stressed")]
analysis_formula = 
  response.c ~  Sex:Drug:Stress:reminder + (1|session) + (Drug:Stress:reminder| participant) + (1|item_name)
if (file.exists(f)) {
  load(f)
} else {
  fit = fit_ordinal_model(
    the_data, analysis_formula, f, family = acat())
}

Click to show sample diagnostics

Table 8.7: Summary of model paramters

parameter	mean	sd	q5	q95	rhat	ess_bulk	ess_tail
Threshold[1]	1.1	0.54	0.3	2.0	1.007	631	753
Threshold[2]	0.3	0.55	-0.6	1.2	1.010	648	918
Threshold[3]	0.1	0.56	-0.8	1.0	1.005	660	956
Threshold[4]	0.3	0.56	-0.7	1.2	1.008	674	1084
Threshold[5]	0.2	0.58	-0.8	1.1	1.006	711	1253
Threshold[6]	-0.1	0.58	-1.0	0.9	1.005	697	1070
Threshold[7]	0.2	0.58	-0.7	1.2	1.006	724	1178
Threshold[8]	0.5	0.59	-0.5	1.5	1.004	768	1340
Threshold[9]	-0.9	0.58	-1.9	0.0	1.007	707	1186
Threshold[10]	0.9	0.57	-0.1	1.8	1.006	698	1089
Threshold[11]	0.3	0.58	-0.6	1.3	1.006	747	898
Threshold[12]	0.0	0.59	-0.9	1.0	1.007	701	1265
Threshold[13]	0.4	0.58	-0.6	1.3	1.004	699	888
Threshold[14]	0.3	0.58	-0.6	1.3	1.006	721	1061
Threshold[15]	1.0	0.60	0.0	2.0	1.007	726	1379
Threshold[16]	0.0	0.59	-1.0	1.0	1.006	753	1050
Threshold[17]	0.5	0.58	-0.5	1.4	1.005	709	998
Threshold[18]	1.3	0.61	0.3	2.3	1.003	757	1502
Threshold[19]	-0.6	0.61	-1.6	0.4	1.006	764	1183
b_Sexman:Drugoxycodone:Stresscontrol:reminderpre	-0.1	0.48	-0.8	0.7	1.010	615	1092
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpre	-0.1	0.47	-0.8	0.8	1.010	617	987
b_Sexman:Drugplacebo:Stresscontrol:reminderpre	0.0	0.48	-0.8	0.8	1.010	615	1067
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpre	-0.1	0.47	-0.8	0.7	1.010	616	977
b_Sexman:Drugoxycodone:Stressstress:reminderpre	0.0	0.48	-0.8	0.8	1.010	614	1128
b_Sexwoman:Drugoxycodone:Stressstress:reminderpre	0.0	0.47	-0.7	0.8	1.010	616	979
b_Sexman:Drugplacebo:Stressstress:reminderpre	0.0	0.48	-0.7	0.8	1.009	613	1077
b_Sexwoman:Drugplacebo:Stressstress:reminderpre	0.1	0.47	-0.7	0.9	1.010	615	1045
b_Sexman:Drugoxycodone:Stresscontrol:reminderpost	0.0	0.47	-0.8	0.8	1.011	619	1096
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpost	-0.1	0.47	-0.8	0.7	1.011	623	1038
b_Sexman:Drugplacebo:Stresscontrol:reminderpost	-0.1	0.48	-0.9	0.7	1.011	616	1095
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpost	-0.1	0.47	-0.9	0.7	1.011	616	985
b_Sexman:Drugoxycodone:Stressstress:reminderpost	0.1	0.48	-0.7	0.9	1.009	621	1128
b_Sexwoman:Drugoxycodone:Stressstress:reminderpost	0.1	0.47	-0.6	0.9	1.011	616	1012
b_Sexman:Drugplacebo:Stressstress:reminderpost	0.1	0.48	-0.7	0.9	1.010	622	1093
b_Sexwoman:Drugplacebo:Stressstress:reminderpost	0.3	0.47	-0.4	1.1	1.011	620	1042
sd_item_name__Intercept	0.5	0.30	0.2	1.2	1.003	1446	1774
sd_participant__Intercept	0.2	0.03	0.2	0.3	1.006	1391	1652
sd_participant__Drugoxycodone:Stresscontrol:reminderpre	0.1	0.04	0.0	0.1	1.004	1010	1664
sd_participant__Drugplacebo:Stresscontrol:reminderpre	0.0	0.03	0.0	0.1	1.001	2222	2015
sd_participant__Drugoxycodone:Stressstress:reminderpre	0.0	0.03	0.0	0.1	1.001	1470	1557
sd_participant__Drugplacebo:Stressstress:reminderpre	0.0	0.03	0.0	0.1	1.000	1458	1865
sd_participant__Drugoxycodone:Stresscontrol:reminderpost	0.1	0.04	0.0	0.1	1.003	1262	1970
sd_participant__Drugplacebo:Stresscontrol:reminderpost	0.0	0.03	0.0	0.1	1.001	2456	2618
sd_participant__Drugoxycodone:Stressstress:reminderpost	0.2	0.04	0.1	0.2	1.001	1644	1711
sd_participant__Drugplacebo:Stressstress:reminderpost	0.2	0.04	0.1	0.2	1.002	1550	1778
sd_session__Intercept	0.0	0.04	0.0	0.1	1.000	1493	1782
model: response.c ~ Sex:Drug:Stress:reminder + (1 | session) + (Drug:Stress:reminder | participant) + (1 | item_name) b_ = fix effect parameters sd_ standard deviation of hierarchical parameters ess = effective sample size

p = plot_buffer_effect(fit,the_data,"stress_emb_selfc")
((p$stress.plot + 
    ylab("stress effect on 'stressed', 'emb.' & 'selfc.' ratings.")) | p$buffer.plot) + 
  plot_annotation(tag_levels = "A") + 
  plot_layout(widths = c(2,3))

Figure 8.10: A: Model-predicted effects of the first stress reminder on ‘stressed’, ‘selfconscious’ and ‘embarrased’ ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Posterior distribution of the interaction of stress effect (pre-post reminder) ratings with drug condition. Dashed vertical lines enclode a region of practical equivalence.

p$buffer.stats %>% 
  kable(caption = "Buffering effect of oxycodone for 'stressed', 'emb.' & 'selfc.' ratings after reminders.") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic()

Table 8.8: Buffering effect of oxycodone for ‘stressed’, ‘emb.’ & ‘selfc.’ ratings after reminders.

Sex	stress	stress in placebo	stress in oxycodone	buffer
men	5 [0, 10]; 0.97975	6 [-1, 13]; 0.9675	3 [-4, 10]; 0.83775	-3 [-13, 7]; 0.28525
women	15 [11, 20]; >.999	21 [14, 28]; >.999	10 [3, 16]; 0.9965	-12 [-22, -2]; 0.00925
total	10 [7, 13]; >.999	14 [9, 18]; >.999	6 [2, 11]; 0.9965	-7 [-14, -0]; 0.023
Sex difference	11 [4, 17]; 0.99875	15 [5, 25]; 0.998	6 [-4, 16]; 0.88525	-9 [-23, 5]; 0.098

8.6 Buffering of negative mood

f = "brmsP1/buffer/negaffect.Rdata"
the_data = my_data.Q_buffer[item_name %in% negaffect_items]
if (file.exists(f)) {
  load(f)
} else {
  fit = fit_ordinal_model(
    the_data, analysis_formula, f, family = acat())
}

Click to show sample diagnostics

Table 8.9: Summary of model paramters

parameter	mean	sd	q5	q95	rhat	ess_bulk	ess_tail
Threshold[1]	1.0	0.60	0.0	2.0	1.002	851	1757
Threshold[2]	0.6	0.61	-0.4	1.6	1.001	864	1748
Threshold[3]	0.9	0.62	-0.1	1.9	1.001	874	1954
Threshold[4]	0.6	0.63	-0.4	1.7	1.002	929	1816
Threshold[5]	0.9	0.64	-0.2	1.9	1.001	922	1916
Threshold[6]	0.8	0.65	-0.3	1.8	1.002	966	1896
Threshold[7]	0.6	0.65	-0.5	1.6	1.002	996	1984
Threshold[8]	1.3	0.67	0.2	2.4	1.001	1052	2125
Threshold[9]	0.4	0.68	-0.8	1.5	1.001	1048	2193
Threshold[10]	1.0	0.67	-0.1	2.1	1.001	1018	2233
Threshold[11]	1.2	0.69	0.0	2.3	1.001	1097	2187
Threshold[12]	0.7	0.69	-0.4	1.9	1.001	1118	2420
Threshold[13]	1.5	0.71	0.3	2.6	1.002	1161	2360
Threshold[14]	1.2	0.75	0.0	2.5	1.001	1277	2810
Threshold[15]	2.2	0.88	0.8	3.7	1.001	1609	3250
Threshold[16]	0.1	0.87	-1.4	1.5	1.001	1632	3193
Threshold[17]	1.9	0.86	0.5	3.3	1.001	1620	3411
Threshold[18]	1.9	1.02	0.3	3.6	1.002	1860	3197
Threshold[19]	-0.1	0.95	-1.7	1.4	1.001	2040	4168
b_Sexman:Drugoxycodone:Stresscontrol:reminderpre	-0.4	0.54	-1.2	0.5	1.003	892	2056
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpre	0.1	0.52	-0.7	1.0	1.003	878	1978
b_Sexman:Drugplacebo:Stresscontrol:reminderpre	-0.4	0.53	-1.2	0.5	1.004	867	1914
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpre	-0.2	0.53	-1.0	0.7	1.004	845	1877
b_Sexman:Drugoxycodone:Stressstress:reminderpre	-0.1	0.52	-1.0	0.7	1.003	861	1862
b_Sexwoman:Drugoxycodone:Stressstress:reminderpre	0.4	0.52	-0.4	1.3	1.003	842	1870
b_Sexman:Drugplacebo:Stressstress:reminderpre	0.0	0.52	-0.9	0.8	1.003	866	2070
b_Sexwoman:Drugplacebo:Stressstress:reminderpre	0.2	0.52	-0.6	1.1	1.004	845	1884
b_Sexman:Drugoxycodone:Stresscontrol:reminderpost	-0.5	0.55	-1.4	0.4	1.003	941	2108
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpost	0.0	0.54	-0.9	0.9	1.003	909	2029
b_Sexman:Drugplacebo:Stresscontrol:reminderpost	-0.6	0.54	-1.4	0.3	1.003	884	1853
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpost	-0.3	0.53	-1.1	0.6	1.004	850	1883
b_Sexman:Drugoxycodone:Stressstress:reminderpost	0.0	0.52	-0.8	0.9	1.004	841	1906
b_Sexwoman:Drugoxycodone:Stressstress:reminderpost	0.6	0.52	-0.2	1.5	1.003	849	1748
b_Sexman:Drugplacebo:Stressstress:reminderpost	0.2	0.52	-0.6	1.1	1.003	848	1913
b_Sexwoman:Drugplacebo:Stressstress:reminderpost	0.6	0.52	-0.2	1.5	1.004	840	1911
sd_item_name__Intercept	0.6	0.38	0.2	1.3	1.003	1492	2131
sd_participant__Intercept	0.7	0.11	0.5	0.8	1.001	2427	4003
sd_participant__Drugoxycodone:Stresscontrol:reminderpre	0.6	0.14	0.3	0.8	1.005	1049	1800
sd_participant__Drugplacebo:Stresscontrol:reminderpre	0.3	0.19	0.0	0.7	1.005	403	1241
sd_participant__Drugoxycodone:Stressstress:reminderpre	0.4	0.12	0.2	0.6	1.003	526	648
sd_participant__Drugplacebo:Stressstress:reminderpre	0.4	0.14	0.2	0.6	1.005	659	753
sd_participant__Drugoxycodone:Stresscontrol:reminderpost	0.7	0.17	0.4	1.0	1.004	1263	2195
sd_participant__Drugplacebo:Stresscontrol:reminderpost	0.4	0.22	0.1	0.8	1.003	434	1152
sd_participant__Drugoxycodone:Stressstress:reminderpost	0.4	0.11	0.2	0.6	1.003	575	1427
sd_participant__Drugplacebo:Stressstress:reminderpost	0.4	0.13	0.2	0.6	1.004	522	662
sd_session__Intercept	0.2	0.23	0.1	0.7	1.001	2363	2642
model: response.c ~ Sex:Drug:Stress:reminder + (1 | session) + (Drug:Stress:reminder | participant) + (1 | item_name) b_ = fix effect parameters sd_ standard deviation of hierarchical parameters ess = effective sample size

p = plot_buffer_effect(fit,the_data,"negaffect")
((p$stress.plot + 
    ylab("stress effect on negative mood ratings")) | p$buffer.plot) + 
  plot_annotation(tag_levels = "A") + 
  plot_layout(widths = c(2,3))

Figure 8.11: A: Model-predicted effects of the first stress reminder on ‘negative mood’, ‘selfconscious’ and ‘embarrased’ ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Posterior distribution of the interaction of stress effect (pre-post reminder) ratings with drug condition. Dashed vertical lines enclode a region of practical equivalence.

p$buffer.stats %>% 
  kable(caption = "Buffering effect of oxycodone for negative mood ratings after reminders.") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic()

Table 8.10: Buffering effect of oxycodone for negative mood ratings after reminders.

Sex	stress	stress in placebo	stress in oxycodone	buffer
men	3 [1, 6]; >.999	5 [2, 8]; >.999	2 [-1, 5]; 0.94825	-2 [-6, 2]; 0.14825
women	8 [6, 11]; >.999	11 [7, 14]; >.999	6 [2, 10]; 0.99875	-5 [-10, 1]; 0.04575
total	6 [4, 8]; >.999	8 [5, 10]; >.999	4 [2, 7]; >.999	-3 [-7, -0]; 0.0235
Sex difference	5 [2, 8]; 0.99825	6 [1, 11]; 0.994625	4 [-1, 9]; 0.93075	-2 [-9, 4]; 0.242125

8.7 Buffering of positive mood

f = "brmsP1/buffer/posaffect.Rdata"
the_data = my_data.Q_buffer[item_name %in% posaffect_items]
if (file.exists(f)) {
  load(f)
} else {
  fit = fit_ordinal_model(
    the_data, analysis_formula, f, family = acat())
}

Click to show sample diagnostics

Table 8.11: Summary of model paramters

parameter	mean	sd	q5	q95	rhat	ess_bulk	ess_tail
Threshold[1]	0.2	0.59	-0.8	1.2	1.007	642	1348
Threshold[2]	-1.5	0.58	-2.5	-0.6	1.007	622	1242
Threshold[3]	-0.9	0.56	-1.8	0.0	1.007	570	1239
Threshold[4]	-0.5	0.55	-1.4	0.4	1.007	572	1077
Threshold[5]	-0.5	0.56	-1.4	0.4	1.008	567	1188
Threshold[6]	-0.5	0.55	-1.4	0.4	1.008	584	1126
Threshold[7]	-0.8	0.55	-1.7	0.1	1.008	541	1078
Threshold[8]	-0.2	0.54	-1.1	0.7	1.008	548	1089
Threshold[9]	-1.2	0.53	-2.0	-0.3	1.008	536	1015
Threshold[10]	0.2	0.53	-0.7	1.1	1.009	514	1047
Threshold[11]	-0.2	0.53	-1.1	0.7	1.008	535	1038
Threshold[12]	-0.2	0.53	-1.1	0.7	1.009	519	1020
Threshold[13]	-0.1	0.53	-1.0	0.7	1.009	520	975
Threshold[14]	0.1	0.53	-0.7	1.0	1.009	519	1037
Threshold[15]	-0.1	0.53	-0.9	0.8	1.008	518	963
Threshold[16]	0.4	0.53	-0.5	1.3	1.008	512	1041
Threshold[17]	0.2	0.53	-0.6	1.1	1.009	517	1031
Threshold[18]	0.6	0.53	-0.2	1.5	1.008	530	1020
Threshold[19]	-0.5	0.53	-1.3	0.4	1.009	516	1028
b_Sexman:Drugoxycodone:Stresscontrol:reminderpre	0.2	0.52	-0.6	1.1	1.009	473	967
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpre	-0.1	0.52	-0.9	0.8	1.007	472	1013
b_Sexman:Drugplacebo:Stresscontrol:reminderpre	0.1	0.51	-0.7	1.0	1.009	475	978
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpre	0.0	0.51	-0.8	0.9	1.007	475	1034
b_Sexman:Drugoxycodone:Stressstress:reminderpre	0.1	0.52	-0.7	1.0	1.009	477	966
b_Sexwoman:Drugoxycodone:Stressstress:reminderpre	-0.2	0.51	-1.1	0.6	1.008	473	1035
b_Sexman:Drugplacebo:Stressstress:reminderpre	0.1	0.51	-0.7	1.0	1.009	476	1002
b_Sexwoman:Drugplacebo:Stressstress:reminderpre	-0.1	0.51	-1.0	0.8	1.007	474	1034
b_Sexman:Drugoxycodone:Stresscontrol:reminderpost	0.2	0.52	-0.6	1.1	1.009	471	972
b_Sexwoman:Drugoxycodone:Stresscontrol:reminderpost	0.0	0.52	-0.8	0.9	1.007	477	1029
b_Sexman:Drugplacebo:Stresscontrol:reminderpost	0.1	0.51	-0.7	1.0	1.009	473	989
b_Sexwoman:Drugplacebo:Stresscontrol:reminderpost	0.0	0.51	-0.8	0.9	1.008	475	1041
b_Sexman:Drugoxycodone:Stressstress:reminderpost	0.0	0.51	-0.8	0.9	1.009	477	1019
b_Sexwoman:Drugoxycodone:Stressstress:reminderpost	-0.3	0.52	-1.1	0.6	1.008	475	1035
b_Sexman:Drugplacebo:Stressstress:reminderpost	0.0	0.51	-0.8	0.9	1.009	475	979
b_Sexwoman:Drugplacebo:Stressstress:reminderpost	-0.3	0.51	-1.1	0.6	1.007	475	1044
sd_item_name__Intercept	0.3	0.14	0.1	0.5	1.001	1752	2514
sd_participant__Intercept	0.3	0.04	0.2	0.4	1.003	901	3238
sd_participant__Drugoxycodone:Stresscontrol:reminderpre	0.3	0.05	0.3	0.4	1.005	1005	1322
sd_participant__Drugplacebo:Stresscontrol:reminderpre	0.1	0.06	0.0	0.2	1.028	170	884
sd_participant__Drugoxycodone:Stressstress:reminderpre	0.3	0.05	0.2	0.4	1.010	463	850
sd_participant__Drugplacebo:Stressstress:reminderpre	0.1	0.05	0.0	0.2	1.016	262	1300
sd_participant__Drugoxycodone:Stresscontrol:reminderpost	0.3	0.05	0.2	0.4	1.004	1133	2156
sd_participant__Drugplacebo:Stresscontrol:reminderpost	0.1	0.06	0.0	0.2	1.028	181	1039
sd_participant__Drugoxycodone:Stressstress:reminderpost	0.3	0.05	0.2	0.4	1.012	362	836
sd_participant__Drugplacebo:Stressstress:reminderpost	0.2	0.05	0.1	0.3	1.022	216	1212
sd_session__Intercept	0.0	0.05	0.0	0.1	1.004	2208	3628
model: response.c ~ Sex:Drug:Stress:reminder + (1 | session) + (Drug:Stress:reminder | participant) + (1 | item_name) b_ = fix effect parameters sd_ standard deviation of hierarchical parameters ess = effective sample size

p = plot_buffer_effect(fit,the_data,"posaffect")
((p$stress.plot + 
    ylab("stress effect on negative mood ratings")) | p$buffer.plot) + 
  plot_annotation(tag_levels = "A") + 
  plot_layout(widths = c(2,3))

Figure 8.12: A: Model-predicted effects of the first stress reminder on ‘negative mood’, ‘selfconscious’ and ‘embarrased’ ratings by state Drug and Sex. Points and lines are means and 66 & 95% Credible intervals. B: Posterior distribution of the interaction of stress effect (pre-post reminder) ratings with drug condition. Dashed vertical lines enclode a region of practical equivalence.

p$buffer.stats %>% 
  kable(caption = "Buffering effect of oxycodone for negative mood ratings after reminders.") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic()

Table 8.12: Buffering effect of oxycodone for negative mood ratings after reminders.

Sex	stress	stress in placebo	stress in oxycodone	buffer
men	-4 [-7, -1]; 0.0075	-6 [-10, -1]; 0.00775	-2 [-7, 2]; 0.163	3 [-3, 10]; 0.847375
women	-8 [-12, -5]; <.001	-10 [-15, -5]; <.001	-7 [-12, -2]; 0.001875	3 [-5, 9]; 0.761375
total	-6 [-9, -4]; <.001	-8 [-11, -4]; <.001	-5 [-8, -1]; 0.003	3 [-2, 8]; 0.883375
Sex difference	-5 [-9, 0]; 0.0325	-4 [-11, 3]; 0.119625	-5 [-11, 2]; 0.0775	-1 [-10, 8]; 0.437125

9 Drug effects in Q5 and Q6

We start again with a look at the data

my_data.Q.drug56 =
  my_data.Q[,.(participant, session, stage, response, state.condition, Sex, Drug, item_name)] %>% 
  .[, session := factor(session)]

stage.names.b = gsub("\\\n"," ", stage.names)
my_data.Q.drug56 %>% 
  .[item_name %in% c("drug_like","drug_again","drug_dislike","high")] %>% 
  .[,.(m = mean(response), se = sd(response)/sqrt(length(unique(participant)))), 
    by = .(stage,Sex,state.condition, Drug, item_name)] %>% 
  .[, `:=`(lower = m-2*se, upper = m+2*se)] %>% 
  ggplot(aes(x = stage, y = m, color = Drug, lty = Sex, fill = Drug)) + 
  geom_line() + 
  facet_grid(item_name~state.condition, scale = "free") + 
  geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .2, color = NA) + 
  scale_col_drug + 
  scale_x_continuous(breaks = 1:11, labels = stage.names.b) + 
  theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust=1))

9.1 Drug liking

Here is a plot of drug liking, measured as the difference between the oxycodone and placebo conditions, dependent on sex.

my_data.Q.drug56 = 
my_data.Q.drug56[stage %in% 5:6]
the_data = my_data.Q.drug56[item_name == "drug_like"]
the_data %>% plot_drug_lines()

Figure 9.1: Drug liking stratified by stress condition and sex. Drug liking is measured as the slope of the lines from placebo to ocxcodone. Each line represents one participant.

We can again summarize the data by plotting the distribution of the difference in liking between drug conditions.

p = plot_drug(the_data)
the_data = the_data %>% coarsen(var = "response", bins = 20, range = c(0,100))
p

Figure 9.2: Distribution of drug liking effects, stratified by stress condition (in panles) and sex (in colors).

raw_contrasts.drug = p$data[, .(participant, Sex, drug_effect, state.condition)][, Outcome := "drug_like"]

Here is the full factorial model we use to estimate the effect of stress:

drug_formula = response.c ~ state.condition:Sex:Drug + (1|session) + (1|participant/stage) + (0+state.condition:Drug|participant)

We estimate the model

fn = "brmsP1/drug/like.Rdata"
if (file.exists(fn)) {
  load(fn)
} else {
  fit = fit_ordinal_model(the_data, drug_formula, fn, family = acat())
}
#pp_check(fit, ndraws = 2000, type = "bars")

Click to show sample diagnostics

Table 9.1: Summary of model paramters

parameter	mean	sd	q5	q95	rhat	ess_bulk	ess_tail
Threshold[1]	0.2	0.97	-1.5	1.8	1.002	772	1375
Threshold[2]	0.3	0.97	-1.3	1.9	1.001	894	1344
Threshold[3]	-0.5	0.99	-2.1	1.1	1.001	878	1497
Threshold[4]	-0.2	0.98	-1.8	1.5	1.002	882	1656
Threshold[5]	0.2	0.98	-1.5	1.8	1.001	887	1451
Threshold[6]	0.3	0.99	-1.3	1.9	1.001	946	1600
Threshold[7]	1.1	1.05	-0.6	2.8	1.001	1089	1765
Threshold[8]	0.0	1.07	-1.8	1.8	1.002	1058	1824
Threshold[9]	0.3	1.04	-1.4	2.0	1.001	994	1655
Threshold[10]	1.1	1.06	-0.7	2.9	1.001	1042	1970
Threshold[11]	-0.2	1.03	-1.9	1.5	1.001	1056	1724
Threshold[12]	1.6	1.05	-0.1	3.3	1.001	1046	1849
Threshold[13]	0.2	1.04	-1.5	1.9	1.001	1089	1884
Threshold[14]	0.9	1.00	-0.8	2.6	1.002	960	1625
Threshold[15]	1.7	1.00	0.1	3.3	1.001	1034	1733
Threshold[16]	1.6	1.03	-0.1	3.3	1.002	967	1895
Threshold[17]	2.6	1.12	0.8	4.5	1.002	1267	2048
Threshold[18]	1.2	1.14	-0.7	3.1	1.001	1317	2535
Threshold[19]	1.0	1.02	-0.7	2.7	1.002	1042	1669
b_state.conditioncontrol:Sexman:Drugplacebo	-1.6	0.89	-3.1	-0.2	1.002	1400	2190
b_state.conditionstress:Sexman:Drugplacebo	-1.0	0.86	-2.5	0.3	1.000	1448	2054
b_state.conditioncontrol:Sexwoman:Drugplacebo	-2.0	0.92	-3.6	-0.6	1.005	1350	2414
b_state.conditionstress:Sexwoman:Drugplacebo	-1.8	0.94	-3.4	-0.3	1.002	1695	2252
b_state.conditioncontrol:Sexman:Drugoxycodone	2.0	0.78	0.8	3.3	1.000	1165	1995
b_state.conditionstress:Sexman:Drugoxycodone	1.8	0.79	0.5	3.1	1.000	1249	1781
b_state.conditioncontrol:Sexwoman:Drugoxycodone	1.3	0.76	0.0	2.6	1.001	1153	1611
b_state.conditionstress:Sexwoman:Drugoxycodone	1.4	0.79	0.1	2.7	1.001	1172	1970
sd_participant__Intercept	0.9	0.28	0.4	1.3	1.010	366	433
sd_participant__state.conditioncontrol:Drugplacebo	1.4	0.53	0.6	2.4	1.003	943	1029
sd_participant__state.conditionstress:Drugplacebo	1.7	0.49	1.0	2.6	1.006	982	1644
sd_participant__state.conditioncontrol:Drugoxycodone	0.8	0.35	0.2	1.4	1.009	323	542
sd_participant__state.conditionstress:Drugoxycodone	1.4	0.34	0.9	2.0	1.005	533	1029
sd_participant:stage__Intercept	0.1	0.07	0.0	0.2	1.002	1558	2235
sd_session__Intercept	1.8	0.87	0.5	3.3	1.003	642	712
model: response.c ~ state.condition:Sex:Drug + (1 | session) + (1 | participant/stage) + (0 + state.condition:Drug | participant) b_ = fix effect parameters sd_ standard deviation of hierarchical parameters ess = effective sample size

For the following plots we calculate the effect of stress on drug liking as liking_in_oxycodone - liking_in_placebo, so that positive numbers mean that liking is greater in the stress condition.

Here are the estimated responses in the “cells” of the design, stratified by sex:

post.contrasts = get_drug_contrasts(fit, the_data)
attr(post.contrasts,"cell_stats") %>% 
  kable(caption = "Estimated responses for 'drug liking' by sex, state condition, drug") %>% 
  kable_styling(full_width = TRUE) %>% 
  kable_classic()

Table 9.2: Estimated responses for ‘drug liking’ by sex, state condition, drug

Sex	state.condition	placebo	oxycodone
men	control	2 (1, 3)	44 (41, 48)
men	stress	7 (5, 9)	41 (38, 43)
women	control	2 (1, 3)	27 (24, 30)
women	stress	3 (2, 4)	34 (31, 36)

post.contrasts = get_drug_contrasts(fit, the_data)
drug_plotter(post.contrasts)

Figure 9.3: Drug liking. Left: Drug liking by sex. Top right: Drug liking by stress condition. Note that drug liking is the difference in liking between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in drug liking, stress effect on drug liking, sex difference in drug liking in the stress condition, sex difference in drug liking in the control condition (for stress). Bottom right: Difference between the sex differences in liking between the control and stress condition.

drug_effects.stats = get_drug_effects.stats(post.contrasts, "drug_like")
drug_effects.stats.b = get_drug_effects.stats.b(post.contrasts, "drug_like")

make_tbl(post.contrasts, 
         caption = "Main and interaction effects of Drug on 'drug liking'.") %>%  
  add_footnote(c("Single words identify main effects in different conditions or sexes.\n '-' identifies contasts. 'in' signifies that an effect is reported for a subgroub or condition.","Numbers show drug effects."), notation="alphabet")

Table 9.3: Main and interaction effects of Drug on ‘drug liking’.

Effecta	mean (CrI)b
women in control	25 [22, 28]; >.999
men in control	42 [39, 46]; >.999
women in stress	31 [28, 34]; >.999
men in stress	34 [31, 37]; >.999
women	28 [26, 30]; >.999
men	38 [36, 40]; >.999
stress	32 [30, 34]; >.999
control	34 [31, 36]; >.999
women - men	-10 [-14, -7]; <.001
stress - ctrl	-1 [-4, 2]; 0.2225
stress - ctrl in men	-9 [-13, -4]; <.001
stress - ctrl in women	6 [2, 11]; 0.99725
women - men in ctrl	-18 [-23, -12]; <.001
women - men in stress	-3 [-7, 2]; 0.1165
three_way	-15 [-22, -8]; <.001
a Single words identify main effects in different conditions or sexes. ‘-’ identifies contasts. ‘in’ signifies that an effect is reported for a subgroub or condition.
b Numbers show drug effects.

9.2 Drug Again

Figure 9.4: Raw data for the effect of drug.

Table 9.4: Estimated responses for drug_again by sex, state condition, drug

Sex	state.condition	placebo	oxycodone
men	control	3 (2, 5)	31 (26, 35)
men	stress	10 (7, 13)	31 (27, 36)
women	control	6 (4, 9)	27 (22, 33)
women	stress	9 (6, 13)	28 (23, 34)

Figure 9.5: Posterior distributions for effect of drug.

The figure shows Left: drug again by sex. Top right: drug again by stress condition. Note that drug again is the difference in drug again between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in drug again, stress effect on drug again, sex difference in drug again in the stress condition, sex difference in drug again in the control condition (for stress). Bottom right: Difference between the sex differences in drug again between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘drug_again’.

Effecta	mean (CrI)b
women in control	21 [15, 27]; >.999
men in control	27 [22, 32]; >.999
women in stress	19 [12, 26]; >.999
men in stress	22 [16, 27]; >.999
women	20 [16, 25]; >.999
men	24 [21, 28]; >.999
stress	20 [16, 25]; >.999
control	24 [20, 28]; >.999
women - men	-4 [-10, 1]; 0.062
stress - ctrl	-4 [-9, 2]; 0.0995
stress - ctrl in men	-6 [-13, 1]; 0.053
stress - ctrl in women	-2 [-11, 8]; 0.368
women - men in ctrl	-6 [-14, 1]; 0.0505
women - men in stress	-2 [-11, 6]; 0.304
three_way	-4 [-16, 7]; 0.2535

9.3 Drug Dislike

Figure 9.6: Raw data for the effect of drug.

Table 9.4: Estimated responses for drug_dislike by sex, state condition, drug

Sex	state.condition	placebo	oxycodone
men	control	5 (3, 8)	11 (7, 15)
men	stress	8 (5, 12)	16 (11, 21)
women	control	3 (2, 5)	24 (19, 29)
women	stress	5 (3, 8)	32 (26, 37)

Figure 9.7: Posterior distributions for effect of drug.

The figure shows Left: drug dislike by sex. Top right: drug dislike by stress condition. Note that drug dislike is the difference in drug dislike between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in drug dislike, stress effect on drug dislike, sex difference in drug dislike in the stress condition, sex difference in drug dislike in the control condition (for stress). Bottom right: Difference between the sex differences in drug dislike between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘drug_dislike’.

Effecta	mean (CrI)b
women in control	21 [16, 26]; >.999
men in control	6 [1, 11]; 0.99475
women in stress	26 [20, 33]; >.999
men in stress	8 [2, 14]; 0.99525
women	24 [20, 28]; >.999
men	7 [3, 10]; >.999
stress	17 [13, 21]; >.999
control	13 [10, 17]; >.999
women - men	17 [11, 22]; >.999
stress - ctrl	4 [-2, 9]; 0.91975
stress - ctrl in men	2 [-6, 10]; 0.705
stress - ctrl in women	6 [-3, 14]; 0.9145
women - men in ctrl	15 [8, 22]; >.999
women - men in stress	19 [10, 27]; >.999
three_way	-4 [-15, 8]; 0.25625

9.4 High

Figure 9.8: Raw data for the effect of drug.

Table 9.4: Estimated responses for high by sex, state condition, drug

Sex	state.condition	placebo	oxycodone
men	control	3 (2, 5)	52 (47, 56)
men	stress	10 (7, 12)	51 (46, 55)
women	control	4 (3, 6)	59 (55, 64)
women	stress	6 (4, 9)	61 (56, 65)

Figure 9.9: Posterior distributions for effect of drug.

The figure shows Left: high by sex. Top right: high by stress condition. Note that high is the difference in high between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in high, stress effect on high, sex difference in high in the stress condition, sex difference in high in the control condition (for stress). Bottom right: Difference between the sex differences in high between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘high’.

Effecta	mean (CrI)b
women in control	55 [50, 60]; >.999
men in control	49 [44, 53]; >.999
women in stress	54 [50, 59]; >.999
men in stress	41 [35, 46]; >.999
women	55 [51, 58]; >.999
men	45 [41, 48]; >.999
stress	48 [44, 51]; >.999
control	52 [48, 55]; >.999
women - men	10 [5, 15]; >.999
stress - ctrl	-4 [-9, 1]; 0.0545
stress - ctrl in men	-8 [-15, -1]; 0.01825
stress - ctrl in women	-0 [-7, 6]; 0.44225
women - men in ctrl	6 [-1, 13]; 0.96575
women - men in stress	13 [6, 21]; >.999
three_way	-7 [-17, 2]; 0.067

9.5 Dizzy

Figure 9.10: Raw data for the effect of drug.

Table 9.4: Estimated responses for dizzy by sex, state condition, drug

Sex	state.condition	placebo	oxycodone
men	control	5 (3, 7)	40 (35, 46)
men	stress	9 (6, 12)	33 (28, 38)
women	control	8 (5, 11)	63 (58, 68)
women	stress	11 (8, 14)	60 (54, 65)

Figure 9.11: Posterior distributions for effect of drug.

The figure shows Left: dizzy by sex. Top right: dizzy by stress condition. Note that dizzy is the difference in dizzy between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in dizzy, stress effect on dizzy, sex difference in dizzy in the stress condition, sex difference in dizzy in the control condition (for stress). Bottom right: Difference between the sex differences in dizzy between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘dizzy’.

Effecta	mean (CrI)b
women in control	55 [49, 61]; >.999
men in control	35 [30, 41]; >.999
women in stress	49 [43, 55]; >.999
men in stress	24 [18, 30]; >.999
women	52 [48, 56]; >.999
men	30 [26, 34]; >.999
stress	37 [32, 41]; >.999
control	45 [41, 49]; >.999
women - men	22 [16, 28]; >.999
stress - ctrl	-9 [-15, -3]; 0.0025
stress - ctrl in men	-11 [-20, -3]; 0.006
stress - ctrl in women	-6 [-14, 3]; 0.082
women - men in ctrl	20 [11, 28]; >.999
women - men in stress	25 [16, 33]; >.999
three_way	-5 [-16, 7]; 0.20875

9.6 Dry Mouth

Figure 9.12: Raw data for the effect of drug.

Table 9.4: Estimated responses for dry_mouth by sex, state condition, drug

Sex	state.condition	placebo	oxycodone
men	control	7 (5, 10)	9 (6, 11)
men	stress	14 (11, 18)	16 (13, 19)
women	control	15 (12, 18)	14 (11, 17)
women	stress	20 (16, 23)	24 (20, 27)

Figure 9.13: Posterior distributions for effect of drug.

The figure shows Left: dry mouth by sex. Top right: dry mouth by stress condition. Note that dry mouth is the difference in dry mouth between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in dry mouth, stress effect on dry mouth, sex difference in dry mouth in the stress condition, sex difference in dry mouth in the control condition (for stress). Bottom right: Difference between the sex differences in dry mouth between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘dry_mouth’.

Effecta	mean (CrI)b
women in control	-1 [-5, 4]; 0.35
men in control	1 [-2, 4]; 0.75675
women in stress	4 [-1, 9]; 0.95025
men in stress	1 [-3, 6]; 0.73775
women	2 [-2, 5]; 0.849
men	1 [-2, 4]; 0.82425
stress	3 [-1, 6]; 0.947
control	0 [-3, 3]; 0.55075
women - men	0 [-4, 5]; 0.5595
stress - ctrl	3 [-2, 7]; 0.88225
stress - ctrl in men	0 [-5, 6]; 0.54675
stress - ctrl in women	5 [-2, 11]; 0.93125
women - men in ctrl	-2 [-7, 3]; 0.2425
women - men in stress	3 [-4, 9]; 0.783
three_way	-5 [-13, 4]; 0.14425

9.7 Blunted

Figure 9.14: Raw data for the effect of drug.

Table 9.4: Estimated responses for blunted by sex, state condition, drug

Sex	state.condition	placebo	oxycodone
men	control	9 [7, 12];	45 [40, 50];
men	stress	14 [10, 17];	35 [30, 40];
women	control	15 [11, 18];	57 [52, 62];
women	stress	14 [11, 17];	53 [48, 58];

Figure 9.15: Posterior distributions for effect of drug.

The figure shows Left: blunted by sex. Top right: blunted by stress condition. Note that blunted is the difference in blunted between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in blunted, stress effect on blunted, sex difference in blunted in the stress condition, sex difference in blunted in the control condition (for stress). Bottom right: Difference between the sex differences in blunted between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘blunted’.

Effecta	mean (CrI)b
women in control	42 [36, 49]; >.999
men in control	36 [30, 42]; >.999
women in stress	39 [33, 45]; >.999
men in stress	22 [16, 27]; >.999
women	41 [37, 45]; >.999
men	29 [25, 33]; >.999
stress	30 [26, 35]; >.999
control	39 [35, 43]; >.999
women - men	12 [6, 18]; >.999
stress - ctrl	-9 [-15, -3]; 0.00325
stress - ctrl in men	-15 [-23, -6]; <.001
stress - ctrl in women	-3 [-12, 6]; 0.2495
women - men in ctrl	6 [-2, 15]; 0.92675
women - men in stress	18 [9, 26]; >.999
three_way	-12 [-24, 0]; 0.02675

9.8 Nauseous

Figure 9.16: Raw data for the effect of drug.

Table 9.4: Estimated responses for nauseous by sex, state condition, drug

Sex	state.condition	placebo	oxycodone
men	control	3 (2, 4)	9 (6, 13)
men	stress	5 (3, 7)	6 (4, 9)
women	control	3 (2, 4)	15 (11, 19)
women	stress	4 (3, 5)	13 (10, 17)

Figure 9.17: Posterior distributions for effect of drug.

The figure shows Left: nauseous by sex. Top right: nauseous by stress condition. Note that nauseous is the difference in nauseous between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in nauseous, stress effect on nauseous, sex difference in nauseous in the stress condition, sex difference in nauseous in the control condition (for stress). Bottom right: Difference between the sex differences in nauseous between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘nauseous’.

Effecta	mean (CrI)b
women in control	12 [8, 16]; >.999
men in control	6 [3, 10]; >.999
women in stress	9 [6, 13]; >.999
men in stress	2 [-2, 5]; 0.823
women	11 [8, 14]; >.999
men	4 [2, 6]; >.999
stress	5 [3, 8]; >.999
control	9 [7, 12]; >.999
women - men	7 [3, 11]; >.999
stress - ctrl	-4 [-7, -0]; 0.0225
stress - ctrl in men	-5 [-10, -0]; 0.0215
stress - ctrl in women	-3 [-9, 3]; 0.1635
women - men in ctrl	6 [1, 11]; 0.98575
women - men in stress	8 [3, 13]; 0.998
three_way	-2 [-10, 6]; 0.32225

9.9 Warm Face

Figure 9.18: Raw data for the effect of drug.

Table 9.4: Estimated responses for warm_face by sex, state condition, drug

Sex	state.condition	placebo	oxycodone
men	control	5 (3, 7)	24 (19, 28)
men	stress	6 (4, 9)	21 (17, 25)
women	control	5 (3, 6)	25 (20, 29)
women	stress	11 (8, 14)	35 (30, 40)

Figure 9.19: Posterior distributions for effect of drug.

The figure shows Left: warm face by sex. Top right: warm face by stress condition. Note that warm face is the difference in warm face between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in warm face, stress effect on warm face, sex difference in warm face in the stress condition, sex difference in warm face in the control condition (for stress). Bottom right: Difference between the sex differences in warm face between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘warm_face’.

Effecta	mean (CrI)b
women in control	20 [16, 25]; >.999
men in control	19 [14, 24]; >.999
women in stress	24 [18, 29]; >.999
men in stress	15 [10, 19]; >.999
women	22 [18, 26]; >.999
men	17 [13, 20]; >.999
stress	19 [16, 23]; >.999
control	20 [16, 23]; >.999
women - men	5 [0, 10]; 0.98075
stress - ctrl	-0 [-6, 5]; 0.43375
stress - ctrl in men	-5 [-11, 2]; 0.0845
stress - ctrl in women	4 [-4, 11]; 0.8295
women - men in ctrl	1 [-6, 8]; 0.6125
women - men in stress	9 [2, 17]; 0.994
three_way	-8 [-18, 2]; 0.05175

9.10 Not Myself

Figure 9.20: Raw data for the effect of drug.

Table 9.4: Estimated responses for not_myself by sex, state condition, drug

Sex	state.condition	placebo	oxycodone
men	control	8 (5, 10)	28 (23, 33)
men	stress	10 (8, 14)	31 (27, 36)
women	control	8 (6, 11)	45 (40, 50)
women	stress	10 (7, 14)	43 (38, 49)

Figure 9.21: Posterior distributions for effect of drug.

The figure shows Left: not myself by sex. Top right: not myself by stress condition. Note that not myself is the difference in not myself between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in not myself, stress effect on not myself, sex difference in not myself in the stress condition, sex difference in not myself in the control condition (for stress). Bottom right: Difference between the sex differences in not myself between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘not_myself’.

Effecta	mean (CrI)b
women in control	37 [31, 42]; >.999
men in control	21 [15, 26]; >.999
women in stress	33 [27, 39]; >.999
men in stress	21 [15, 27]; >.999
women	35 [31, 39]; >.999
men	21 [17, 25]; >.999
stress	27 [23, 31]; >.999
control	29 [25, 33]; >.999
women - men	14 [8, 20]; >.999
stress - ctrl	-2 [-7, 4]; 0.28475
stress - ctrl in men	0 [-7, 8]; 0.53875
stress - ctrl in women	-4 [-12, 5]; 0.19475
women - men in ctrl	16 [8, 24]; >.999
women - men in stress	12 [4, 20]; 0.99825
three_way	4 [-7, 15]; 0.75825

9.11 Pos Affect

Figure 9.22: Raw data for the effect of drug.

Table 9.4: Estimated responses for pos_affect by sex, state condition, drug

Sex	state.condition	placebo	oxycodone
men	control	72 (70, 74)	73 (71, 75)
men	stress	68 (66, 70)	69 (67, 71)
women	control	68 (66, 70)	64 (62, 66)
women	stress	57 (55, 59)	53 (51, 55)

Figure 9.23: Posterior distributions for effect of drug.

The figure shows Left: pos affect by sex. Top right: pos affect by stress condition. Note that pos affect is the difference in pos affect between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in pos affect, stress effect on pos affect, sex difference in pos affect in the stress condition, sex difference in pos affect in the control condition (for stress). Bottom right: Difference between the sex differences in pos affect between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘pos_affect’.

Effecta	mean (CrI)b
women in control	-4 [-7, -1]; 0.00175
men in control	1 [-2, 3]; 0.66475
women in stress	-4 [-7, -1]; 0.0055
men in stress	1 [-2, 4]; 0.789
women	-4 [-6, -2]; <.001
men	1 [-1, 3]; 0.81025
stress	-1 [-3, 1]; 0.08075
control	-2 [-4, 0]; 0.03775
women - men	-5 [-8, -2]; <.001
stress - ctrl	0 [-2, 3]; 0.57825
stress - ctrl in men	1 [-3, 4]; 0.60825
stress - ctrl in women	0 [-4, 4]; 0.508
women - men in ctrl	-5 [-8, -1]; 0.0085
women - men in stress	-5 [-9, -1]; 0.00575
three_way	0 [-5, 6]; 0.56825

9.12 Neg Affect

Figure 9.24: Raw data for the effect of drug.

Table 9.4: Estimated responses for neg_affect by sex, state condition, drug

Sex	state.condition	placebo	oxycodone
men	control	4 (3, 5)	5 (4, 6)
men	stress	10 (8, 11)	9 (7, 10)
women	control	5 (4, 6)	9 (8, 10)
women	stress	19 (17, 22)	18 (16, 20)

Figure 9.25: Posterior distributions for effect of drug.

The figure shows Left: neg affect by sex. Top right: neg affect by stress condition. Note that neg affect is the difference in neg affect between the ocycodone and placebo conditions. Bottom right, clockwise from top left: Sex differences in neg affect, stress effect on neg affect, sex difference in neg affect in the stress condition, sex difference in neg affect in the control condition (for stress). Bottom right: Difference between the sex differences in neg affect between the control and stress condition.

Table 9.4: Main and interaction effects of drug on ‘neg_affect’.

Effecta	mean (CrI)b
women in control	4 [2, 6]; >.999
men in control	1 [-0, 2]; 0.962
women in stress	-1 [-4, 2]; 0.1995
men in stress	-1 [-3, 1]; 0.1375
women	1 [-0, 3]; 0.94125
men	0 [-1, 1]; 0.51775
stress	-1 [-3, 1]; 0.1005
control	3 [2, 4]; >.999
women - men	1 [-1, 3]; 0.90625
stress - ctrl	-4 [-6, -2]; <.001
stress - ctrl in men	-2 [-5, 0]; 0.031
stress - ctrl in women	-5 [-8, -2]; <.001
women - men in ctrl	3 [1, 5]; 0.99625
women - men in stress	-0 [-4, 3]; 0.45525
three_way	3 [-1, 7]; 0.93375

9.13 Summary of Drug effects

drug_effects.stats %>% 
  .[Outcome == "high", Outcome := "feeling_high"] %>% 
  .[grepl("high|drug|dizzy|blunted|not|warm", Outcome)] %>% 
  .[,Outcome := stringr::str_to_title(Outcome)] %>% 
  .[, Outcome := gsub("_"," ", Outcome)] %>% 
  kable(caption = "Drug effects on subjective ratings.") %>% 
  kable_styling(full_width = FALSE) %>% 
  kable_classic() 

Table 9.5: Drug effects on subjective ratings.

Outcome	Men	Women	Total	Sex difference
Drug like	38 [36, 40];	28 [26, 30];	33 [31, 34];	-10 [-14, -7]; <.001
Drug again	24 [21, 28];	20 [16, 25];	22 [19, 25];	-4 [-10, 1]; 0.062
Drug dislike	7 [3, 10];	24 [20, 28];	15 [13, 18];	17 [11, 22]; >.999
Feeling high	45 [41, 48];	55 [51, 58];	50 [47, 52];	10 [5, 15]; >.999
Dizzy	30 [26, 34];	52 [48, 56];	41 [38, 44];	22 [16, 28]; >.999
Blunted	29 [25, 33];	41 [37, 45];	35 [32, 38];	12 [6, 18]; >.999
Warm face	17 [13, 20];	22 [18, 26];	19 [17, 22];	5 [0, 10]; 0.98075
Not myself	21 [17, 25];	35 [31, 39];	28 [25, 31];	14 [8, 20]; >.999

clean = function(dt) {
  dt = 
    dt %>%  
    .[grepl("lik|high", Outcome)] %>%  
    .[Outcome == "high", Outcome := "feeling_high"] %>% 
    .[,Outcome := stringr::str_to_title(Outcome)] %>%
    .[,Sex := ifelse(Sex == "men","Men","Women")] %>%
    .[,state.condition := stringr::str_to_title(state.condition)] %>%
    .[, Outcome := gsub("_"," ", Outcome)] %>% 
    .[, Outcome := gsub("like","liking", Outcome)] %>% 
    .[, Outcome := paste("\u0394 ",Outcome)] %>% 
    .[, Outcome := factor(Outcome)]
  o.levels = levels(dt$Outcome)[c(3,2,1)]
  dt[, Outcome := factor(Outcome, levels = o.levels)]
  return(dt)
}

raw_contrasts.drug.p = 
  raw_contrasts.drug %>%
  copy() %>% 
  clean() %>%  
  .[, l := fmean(drug_effect)-3*fsd(drug_effect), by = .(Outcome,Sex)] %>% 
  .[, u := fmean(drug_effect)+3*fsd(drug_effect), by = .(Outcome,Sex)] %>% 
  .[, plot := drug_effect < u & drug_effect > l]

breaks_drug_fun <- function(x) {
  d.count <<- d.count + 1L
  switch(
    d.count,
    seq(0,100,25), c(0),
    seq(-25,100,25),c(0),
    seq(-50,100,50), c(0),
    seq(0,100,25), c(0),
    seq(-25,100,25),c(0),
    seq(-50,100,50),
  )
}

d.count = 0
# summary_plot.drug = 
#   drug_effects.stats.b %>% 
#     clean() %>% 
#     ggplot(aes(x = state.condition, y = drug_effect,
#                fill = Sex, shape = state.condition,
#                group = factor(Sex):factor(state.condition))) + 
#     geom_bar(stat = "identity", position = position_dodge(1)) + 
#     geom_hline(yintercept = 0) +
#     geom_errorbar(aes(ymin = lower, ymax = upper),  linewidth = 1,
#                   width = 1/5, position = position_dodge(1)) + 
#     geom_quasirandom(data = raw_contrasts.drug.p[plot == TRUE],
#                      dodge.width = 1,
#                      color = adjustcolor("black",alpha = .5),
#                      aes(shape = state.condition)) +
#     scale_shape_manual(values = c(1,2)) + 
#     facet_grid(Outcome ~ Sex, scale = "free_y",switch = "y") + 
#     #scale_y_continuous(breaks = breaks_drug_fun()) + 
#     guides(fill = "none", 
#            shape = guide_legend(nrow = 1, title.position = "left", title = "State induction")) + 
#     theme(legend.position = "top", axis.line.y = element_line(),
#           panel.grid.major.y = element_blank(), panel.spacing = unit(1, "lines"),
#           strip.background = element_blank(), legend.key.size = unit(1,"line"),
#           strip.placement = "outside", legend.title=element_text(size=12)) + 
#     scale_col_sex2 + ylab("") + xlab("")
summary_plot.drug = 
  drug_effects.stats.b %>% 
  clean() %>% 
  ggplot(aes(x = Sex, y = drug_effect,
             fill = Sex, shape = state.condition,
             group = factor(Sex):factor(state.condition))) + 
  geom_bar(stat = "identity", position = position_dodge(1)) + 
  geom_hline(yintercept = 0) +
  geom_errorbar(aes(ymin = lower, ymax = upper),  linewidth = 1,
                width = 1/5, position = position_dodge(1)) + 
  geom_quasirandom(data = raw_contrasts.drug.p[plot == TRUE],
                   dodge.width = 1,
                   color = adjustcolor("black",alpha = .5),
                   aes(shape = state.condition)) +
  scale_shape_manual(values = c(1,2)) + 
  facet_grid(Outcome ~ ., scale = "free_y",switch = "y") + 
  #scale_y_continuous(breaks = breaks_drug_fun()) + 
  guides(fill = "none", 
         shape = guide_legend(nrow = 1, title.position = "left", title = "State induction")) + 
  theme(legend.position = "top", axis.line.y = element_line(),
        panel.grid.major.y = element_blank(), panel.spacing = unit(1, "lines"),
        strip.background = element_blank(), legend.key.size = unit(1,"line"),
        strip.placement = "outside", legend.title=element_text(size=12)) + 
  scale_col_sex2 + ylab("") + xlab("") + gg_no_x_axis

emb_plot = 
  select_data("embarassed") %>% 
    plot_stages(
        my_stages = 2:9,
        by.Sex = TRUE, lty = "Sex",
        shape = "state.condition", 
        by.Drug = TRUE, fill = "Drug") + 
    facet_wrap(~Sex, nrow = 2) + 
    theme(strip.background = element_blank(),
          legend.position = c(.5,.9),
          legend.margin=margin(0,0,0,0),
          legend.spacing = unit(.25,"line"),
          legend.title=element_text(size=12)) + 
    coord_cartesian(ylim = c(0,100), xlim = c(-20,35)) + 
    guides(fill = guide_legend(nrow = 2, title.position = "top", title = "Drug"), 
           shape = "none",
           lty = "none", color = "none") +
    ylab("Embarrassed (VAS)")

summary_plot = 
  (summary_plot.stress | summary_plot.drug | emb_plot) + 
  plot_layout(widths = c(1,2,1.5)) +
  plot_annotation(tag_levels = "A")
s.count = 0
#d.count = 0
ggsave(summary_plot,file = "figures/summary_plot.png", width = 9, height = 7)
s.count = 0
#d.count = 0
ggsave(summary_plot,file = "figures/summary_plot.pdf", width = 9, height = 7)

Overview over effects of stress and drug

---

1. For instance, when examining the stress effect dependent of drug, we have 8 conditions (2 drug conditions time 2 stress conditions times 2 measurements pre and post stress induction). Then we calculate for each of these 8 conditions the response probability for each of the 20 response levels.↩︎